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Abstract 

The  Space  Surveillance  Telescope  (SST)  is  a  Defense  Advanced  Research 
Projects  Agency  (DARPA)  program  designed  to  detect  objects  in  space  like  Near  Earth 
Asteroids  (NEAs)  and  space  debris  in  the  Geosynchronous  Earth  Orbit  (GEO)  belt. 
Binary  hypothesis  tests  (BHTs)  have  historically  been  used  to  facilitate  the  detection  of 
new  objects  in  space.  In  this  dissertation,  a  multi-hypothesis  test  (MHT)  detection 
strategy  is  introduced  to  improve  the  detection  performance  of  the  SST.  In  this  context, 
the  MHT  detennines  if  an  unresolvable  point  source  is  in  the  center,  comer  or  side  of  a 
pixel  in  contrast  to  a  BHT,  which  only  tests  whether  an  object  is  in  the  pixel  or  not.  An 
experiment,  recording  observations  of  a  known  GEO  satellite  as  it  enters  eclipse,  is  used 
to  demonstrate  improved  probability  of  detection  with  the  MHT  by  as  much  as  50%  over 
existing  BHT  methods. 

In  order  to  achieve  optimal  perfonnance  of  the  SST,  alignment  of  the  telescope  is 
conducted  by  retrieving  phase  information  from  defocused  point  sources  to  detennine  the 
telescope’s  aberrations  and  then  the  mirrors  are  moved  for  optical  correction.  A  new 
direct  search  phase  retrieval  technique  for  detennining  the  optical  prescription  of  an 
imaging  system  in  terms  of  Zernike  coefficients  is  described.  The  technique  provides 
coefficient  estimates  without  the  need  to  defocus  point  source  images  to  generate  phase 
diversity  by  using  electric  field  estimates  in  addition  to  intensity  data.  Simulated  point 
source  data  shows  the  new  phase  retrieval  algorithm  avoids  getting  trapped  in  local 
minima  over  a  wide  range  of  random  aberrations.  Experimental  point  source  data  are 
used  to  demonstrate  the  phase  retrieval  effectiveness. 
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ENHANCING  GROUND  BASED  TELESCOPE  PERFORMANCE  WITH  IMAGE 


PROCESSING 

I.  Introduction 

The  Department  of  Defense  recently  fielded  an  f/1  Mersenne-Schmidt  telescope 
called  the  Space  Surveillance  Telescope  (SST),  which  saw  first  light  on  15  February  2011 
[1].  The  SST  has  significantly  advanced  the  United  State’s  ability  to  maintain  space 
situational  awareness  (SSA)  beyond  that  provided  by  the  operationally  employed 
Ground-based  Electro-Optical  Deep  Space  Surveillance  (GEODSS)  telescopes  [2].  The 
main  advantages  of  the  SST  are  that  it  has  a  3.5  m  primary  mirror  and  a  6  deg  wide  field- 
of-view  (FOV)  [3].  In  contrast,  the  GEODSS  f/2.15  Ritchey-Chretien  designed 
telescopes  only  have  aim  primary  mirror  and  a  1 .68  deg  FOV  [2].  The  larger  primary 
mirror  and  increased  FOV  allow  the  SST  to  scan  a  larger  portion  of  the  sky  in  a  shorter 
period  of  time  with  improved  detection  perfonnance  over  GEODSS. 

SSA  is  a  critical  military  mission  and  it  directly  supports  the  US  National  Space 
Policy  to  “(p)reserve  the  Space  Environment... the  United  States  shall  develop,  maintain, 
and  use  space  situational  awareness  infonnation  from  commercial,  civil,  and  national 
security  sources  to  detect,  identify,  and  attribute  actions  in  space  that  are  contrary  to 
responsible  use  and  the  long-tenn  sustainability  of  the  space  environment  [4].”  The  SST 
fills  an  important  niche  in  the  nation’s  space  surveillance  network  by  providing  timely 
and  accurate  updates  to  the  Joint  Space  Operations  Center’s  (JSpOC’s)  space  catalogue 
[5,  6,  7].  Through  synoptic  search  of  deep  space  (i.e.  GEosynchronous  Orbit  (GEO)  and 
Highly  Elliptical  Orbit  (HEO))  on  a  regular  basis,  the  SST  can  detect  and  detennine  the 
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orbits  of  previously  unknown  space  objects.  These  previously  unknown  space  objects  are 
commonly  called  uncorrelated  targets  (UCTs).  Timely  updates  of  these  UCTs  to  the 
JSpOC’s  space  catalogue  support  vital  decision  making  by  US  Strategic  Command’s 
Joint  Functional  Component  Command  for  Space  in  support  of  the  US  National  Space 
Policy. 

Three  major  threats  from  space  to  both  our  satellite  networks  and  the  earth  have 
been  identified.  These  threats  include:  space  debris,  micro-satellites,  and  near  earth 
asteroids  (NEAs)  [8,  9,  10].  The  detection  and  characterization  of  these  threats  can 
provide  the  early  warning  necessary  to  take  any  responsive  actions.  While  the  SST  was 
designed  for  the  military  mission  of  detecting  debris  and  microsatellites  in  deep  space,  it 
is  also  being  used  for  the  detection  of  other  astronomical  objects  like  NEAs.  The  SST’s 
NEA  detection  work  is  being  done  in  partnerships  with  the  US  Naval  Observatory 
(USNO)  and  the  National  Air  and  Space  Administration  (NASA).  However,  the  US  Air 
Force’s  primary  concern  is  the  protection  of  critical  national  space  assets  in  earth  orbit 
from  space  debris  and  micro-satellites  [11]. 

Figure  1  illustrates  the  known  objects  (both  satellites  and  debris)  that  are 
cataloged  by  NASA  and  highlights  the  sheer  number  of  objects  currently  being  tracked. 
The  SSA  functions  that  are  critical  to  avoiding  collisions  of  these  objects  in  space  consist 
of  the  detection  of  UCTs,  accurately  detennining  their  orbits,  and  maintaining  an  up-to- 
date  space  catalogue.  The  SST’s  main  roles  in  SSA  are  the  detection  and  orbit 
detennination  of  UCTs  in  deep  space.  GEO  can  be  distinguished  in  Figure  1  as  the  dense 
ring  of  objects  near  the  equatorial  plane.  In  contrast,  radar  system  like  the  space  fence 
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are  better  suited  for  finding  objects  in  Low  Earth  Orbit  (LEO),  which  is  the  dense  cloud 
of  objects  near  the  earth  in  Figure  1  [12]. 


Figure  1.  A  NASA  produced  image  depicting  the  number  of  satellites  and  debris  tracked  in 
earth’s  orbit.  (Note  that  the  objects  are  not  to  scale)  [12] 


A  3-D  optical  design  layout  and  scale  image  of  SST  is  shown  in  Figure  2.  The 
driving  design  requirements  for  the  telescope  differ  from  more  typical  astronomical 
telescopes.  The  main  difference  is  that  SST  needs  to  be  able  to  scan  deep  space  on  a 
regular  basis  to  detect  and  track  UCTs  versus  maintaining  accurate  orbits  of  known 
objects.  In  order  to  accomplish  that  SSA  mission,  the  Mersenne-Schmidt  design  was 
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selected  for  both  its  wide  FOV  and  compact  design  [13].  The  size  of  the  primary  mirror 
was  driven  by  the  need  to  detect  small  faint  objects  with  relatively  short  integration  times 
to  avoid  streaking  of  the  satellite  image  across  multiple  charged  coupled  device  (CCD) 
pixels.  One  drawback  to  the  design  is  that  the  optical  wavefront  is  curved  in  the  image 
plane  due  to  telescope  aberrations.  To  alleviate  those  optical  aberrations,  a  unique  curved 
detector  array  was  fabricated  and  implemented  in  the  SST  camera  [14]. 


Figure  2.  An  optical  design  schematic  and  scale  picture  of  the  SST. 


1.1  Motivation 

This  dissertation  investigates  how  image  processing  can  improve  the  SST’s 
detection  capabilities  without  requiring  physical  design  changes  to  the  telescope.  The 
two  objectives  that  were  introduced  in  the  prospectus  for  this  work  have  been  met  and  are 
listed  below. 
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Objective  1.  Characterize  the  SST  using  a  novel  aberration  estimator 

Objective  2.  Enhance  the  SST’s  perfonnance  using  detection  and  estimation  theory 

The  SST  was  designed  and  built  by  Defense  Advanced  Research  Projects 
Agency’s  (DARPA)  and  has  completed  its  system  test  and  demonstration  (T&D).  For 
that  phase,  a  set  of  metrics  were  defined  to  evaluate  the  telescope’s  ability  to  perform  the 
SSA  mission.  According  to  the  SST  Phase  I  System  T&D  Procedure,  “(e)ach  system 
perfonnance  metric  (SPM)  is  traceable  to  component  and  subsystem  specifications  of  the 
underlying  technology  development. .  .(t)hey  are  the  core  metrics  that  constitute  the  crux 
of  the  SST’s  performance,  those  that  will  enable  the  efficiencies  of  synoptic  search  [15].” 
Therefore,  improvements  to  these  SPMs  will  provide  more  overall  system  capability.  The 
metrics  are: 

I.  Search  rate:  Q  [deg“/hr]  vs.  mccd  (instrument  magnitude) 

II.  Metric  accuracy:  AO  [arc  sec]  at  nominal  mccd 

III.  Photometric  accuracy:  A  mcc d  at  nominal  mccd 

IV.  Sensitivity:  mccd  vs.  t  (integration  time)  in  various  tracking  modes 

Metric  I  measures  the  search  rate,  Q,  with  different  integration  times  (i.e.  target 
instrument  magnitude  (mccd))  to  quantify  the  performance  traded  off  between  sensitivity 
and  search  rate.  Metric  II  measures  the  telescope’s  metric  accuracy,  AO,  as  a  function  of 
Weed-  Metric  accuracy  represents  error  in  the  estimates  of  an  object’s  azimuth  and 
elevation  in  units  of  arc  seconds.  Metric  III  evaluates  the  variance  of  mccd  and  is 
produced  by  comparing  the  object  of  interest  brightness  (i.e.  digital  counts)  to  the  known 
brightness  of  stars  in  the  FOV.  Metric  IV,  the  system  sensitivity,  is  measured  in  tenns  of 
mccd  and  is  evaluated  as  a  function  of  integration  time,  t. 

The  remote  sensing  research  accomplished  and  discussed  in  this  dissertation  can 
improve  the  SST’s  performance  on  three  of  the  four  SPMs.  In  particular,  the  multi- 
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hypothesis  test  (MHT)  that  is  introduced  in  Chapter  IV  has  new  properties  that  could  be 
used  as  the  kernel  of  an  entirely  new  image  processing  scheme.  This  algorithm  has  the 
potential  to  significantly  outperform  the  baseline  algorithm.  Metric  accuracy  can  be 
improved  with  the  sub-pixel  position  infonnation  inherent  to  the  MHT.  Photometric 
accuracy  should  also  be  improved  because  the  output  of  the  MHT  can  be  used  to  estimate 
photons  in  both  the  pixel  in  which  the  object  is  detected  and  the  neighboring  pixels. 
Finally  and  most  importantly,  the  sensitivity  gains  provided  by  the  MHT  are 
demonstrated  using  experimental  data  from  the  SST. 

1.2  Accomplished  Work 

Improving  the  detection  perfonnance  of  the  SST  through  image  processing 
enhances  the  detection  of  dim  objects  and  the  rejection  of  false  alanns  without  requiring 
significant  changes  to  the  telescope’s  design  or  operations  [16].  The  following  three 
concepts  were  identified  in  the  prospectus  and  investigated  in  the  research  for  improving 
the  detection  sensitivity  of  the  SST. 

1.  Determine  a  better  way  to  retrieve  phase  information  from  the  SST  data. 

2.  Investigate  multi -hypothesis  testing  (MHT)  to  improve  detection  performance. 

3.  Mitigate  star  crossings  from  affecting  the  detection  of  an  object. 

This  document  is  divided  into  six  chapters  that  provide  background,  details  on  the 
work  accomplished,  and  final  conclusions  that  are  intended  to  show  how  these  concepts 
were  investigated.  In  Chapter  II,  the  background  information  provided  includes  a 
literature  review,  an  introduction  to  the  fundamental  image  processing  concepts  used  in 
this  work,  and  the  foundation  for  this  research.  There  are  two  chapters  on  phase  retrieval. 
Chapter  III  describes  a  phase  retrieval  technique  with  a  long  exposure  atmospheric  model 
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useful  for  the  current  SST  shutter  camera,  while  in  Chapter  V,  phase  retrieval  with  a  short 
exposure  model  is  investigated  for  a  future  frame  transfer  camera.  Phase  infonnation  is 
important  because  it  is  used  to  quantify  the  presence  of  optical  aberrations  in  the 
telescope  [17].  More  accurate  infonnation  about  the  phase  should  help  improve  the 
process  of  focus  and  alignment  resulting  in  sharper  images.  Focusing  reduces  the  point 
spread  function  (PSF)  and  increases  the  signal-to-noise  ratio  (SNR)  so  that  the  SST’s 
detection  performance  will  be  improved  using  its  existing  and  future  detection 
algorithms.  The  Cramer-Rao  Lower  Bound  (CRLB)  calculations  demonstrate  that 
improved  aberration  estimates  are  possible  in  certain  coefficient  ranges.  New  algorithms 
are  then  derived  to  accomplish  improved  optical  aberration  estimates  [18,  19]. 

Both  long  and  short  exposure  atmospheric  models  are  used  in  the  different  phase 
retrieval  techniques  investigated  [20],  The  long  exposure  model  is  necessary  for  the 
current  camera  and  physical  shutter  on  the  SST  system.  However,  future  variants  of  the 
SST  or  other  ground-based  three  mirror  telescopes  like  the  Large  Synoptic  Survey 
Telescope  (LSST)  may  benefit  from  phase  retrieval  with  a  short  exposure  atmospheric 
model  [21].  For  the  SST,  if  a  new  frame  transfer  camera  is  procured  it  would  enable  the 
telescope  to  take  short  exposure  images.  Both  long  and  short  exposure  phase  retrieval 
techniques  have  performance  advantages  and  limitations  that  are  discussed  in  Chapter  VI. 

Phase  retrieval  is  also  critical  to  achieving  the  potential  detection  and  false  alarm 
improvements  afforded  by  the  new  detection  algorithm  described  in  Chapter  IV.  The 
new  detection  algorithm  leverages  the  long  exposure  phase  retrieval  outputs  to  form  the 
PSF  estimate  used  in  a  new  multi-hypothesis  test  (MHT)  for  the  detection  of  dim  objects. 
Choosing  the  best  available  detection  algorithm  is  critical  to  maximizing  detection 
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sensitivity  of  the  telescope.  A  thorough  comparison  of  the  current  and  proposed 
detection  algorithms  using  satellite  eclipse  experiment  data  provides  a  side-by-side 
comparison  of  techniques.  Algorithm  perfonnance  assessment  is  based  on  detection 
improvement,  computational  burden  and  implementation  complexity.  The  MHT 
algorithm  mitigates  the  effects  of  aliasing  caused  by  undersampling  of  the  image  and 
employs  a  new  method  of  background  noise  calculations  that  reduces  the  noise 
degradation  associated  with  star  crossings.  The  MHT  significantly  improves  the 
probability  of  detecting  uncorrelated  targets  (UCTs)  over  the  algorithm  currently  used  by 
the  SST.  The  Chapter  VI  conclusions  highlight  the  significant  improvement  in  the  SST’s 
detection  sensitivity  and  the  phase  retrieval  findings  discovered  in  this  research  effort. 
Suggestions  for  future  areas  of  investigation  to  further  enhance  the  perfonnance  of  the 
telescope  are  also  discussed. 
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II.  Background 


The  act  of  finding  space  objects  with  ground  based  telescopes  has  been  performed 
by  astronomers  since  the  17th  century  [22].  The  incorporation  of  CCD  arrays  in  place  of 
film  or  the  human  eye  has  enabled  modem  image  processing  techniques  to  be  used  for 
detection  of  space  objects  [23,  24,  25].  The  SST’s  current  detection  method  is  based  on 
an  algorithm  used  for  the  Lincoln  Near  Earth  Asteroid  Research  (LINEAR)  mission 
conducted  at  the  Experimental  Test  Site  near  Socorro,  NM.  Viggh  et  al.  described  the 
LINEAR  detection  algorithm  using  a  similar  block  diagram  to  the  one  shown  in  Figure  3. 


Figure  3.  The  SST  detection  block  diagram  [25] 


The  SST’s  current  demonstration  software  uses  input  data  that  is  comprised  of  three  to 
five  image  frames  of  the  sky  in  either  a  sidereal  or  a  satellite  track  mode.  Each  frame  is 
exposed  over  a  user  defined  integration  time  period  optimized  to  find  UCTs.  Image 
registration  corrects  for  telescope  pointing  errors  by  using  stars  with  known  coordinates 
in  each  of  the  image  frames.  The  single  frame  point  detection  of  an  object  is  performed 
in  a  two  step  process  of  background  suppression  normalization  followed  by  binary 
quantization,  which  is  described  mathematically  in  Section  2.4.  In  each  frame,  adjacent 
pixels  that  have  objects  above  the  detection  threshold  are  clustered  together  and  classified 
as  a  single  object.  The  algorithm  then  determines  the  centroid  and  extent  of  each 


9 


clustered  object.  The  data  is  then  filtered  using  a  velocity  matched  filter  to  remove  any 
objects  not  moving  at  the  rate  of  the  object(s)  of  interest.  This  is  possible  because  in  rate 
track  mode  the  object  in  Earth’s  orbit  will  be  stationary  and  the  stars  will  be  moving  in 
the  image  between  frames  at  the  sidereal  rate.  The  opposite  is  true  when  the  telescope  is 
in  sidereal  track  mode.  In  addition,  in  order  to  reduce  false  alarms  without  changing  the 
single  frame  detection  threshold,  objects  that  are  not  detected  in  three  to  five  consecutive 
frames  are  also  removed.  Once  the  UTC  is  identified,  the  visual  magnitude  of  that  object 
is  detennined  with  a  plate  model  that  uses  calibration  stars  in  the  FOV  [3,  26].  Finally, 
the  orbit  determination  of  the  UCT  is  determined  by  revisiting  the  object’s  initial 
coordinates  at  a  later  time. 

The  critical  step  in  the  existing  software  described  above  is  the  single  frame  point 
detection  step  because  only  objects  that  exceed  the  detection  threshold  have  the  potential 
for  being  discovered  by  the  detection  algorithm.  Therefore,  improvements  in  single 
frame  point  detection  will  improve  the  overall  perfonnance  of  the  system.  Two  possible 
methods  for  improving  the  single  frame  detection  are  reducing  the  spot  size  to  increase 
the  pixel  SNR  or  to  use  a  matched  filtering  operation  based  on  a  model  of  the  PSF.  The 
phase  retrieval  methods  discussed  in  this  dissertation  have  the  potential  to  improve  the 
understanding  of  the  telescope’s  aberrations  for  focus  and  alignment.  These  methods 
could  be  used  to  decrease  spot  size  and  the  detection  methods  are  proven  to  provide 
enhanced  single-frame  probability  of  detection. 

The  background  provided  in  this  section  is  intended  to  cover  the  fundamental 
concepts  reviewed  in  literature  pertinent  to  the  completed  research.  The  telescope  model, 
Zernike  polynomials  and  atmospheric  models  described  are  well  accepted  in  the  field  of 
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applied  optics  and  form  the  underlying  principles  for  the  new  techniques  developed  in  the 
phase  retrieval  and  detection  chapters. 

2.1  Telescope  Model 

For  both  phase  retrieval  and  the  MHT,  a  wave  model  of  the  telescope  is  used  to 
describe  the  optical  aberrations  and  resulting  PSF  model.  The  aberrations  are 
parameterized  using  Zemike  polynomials  because  they  are  directly  related  to  the 
wavefront  errors  associated  with  the  optical  prescription  of  the  telescope  [17].  In 
addition,  the  Zemike  polynomials  can  be  used  to  model  the  atmosphere  for  short 
exposure  imagery  [27,  28].  The  polynomials  are  defined  on  the  unit  circle  and  form  an 
orthonormal  basis  set  for  polynomial  decomposition  of  wavefront  error  [29]. 

2.1.1  Zemike  Polynomial  for  Defocus 

One  of  the  easiest  ways  to  understand  how  wavefront  errors  are  introduced  into  an 
optical  system  is  by  studying  the  aberration  known  as  defocus.  In  an  ideal  single  lens 
optical  system,  the  rays  of  light  associated  with  a  plane  wave  produced  by  an  object 
infinitely  far  from  the  lens  (effectively  a  point  source)  are  focused  to  a  diffraction-limited 
spot  at  the  focal  point  of  the  lens  as  illustrated  in  Figure  4  [30,  31].  Moving  the  image 
plane  before  or  after  the  focal  plane  causes  defocus  to  occur,  this  introduces  a  wavefront 
error  in  the  pupil  plane  that  can  be  quantified  using  the  Zernike  polynomial  for  defocus. 
The  farther  the  image  plane  is  from  the  focal  plane,  the  more  the  Zernike  polynomial  for 
defocus  has  to  be  scaled  using  the  Zernike  coefficient  for  defocus  to  accurately  model  the 
additional  wavefront  error  of  the  optical  system. 
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Plane  Wave 


Figure  4.  Illustration  depicting  how  moving  the  image  plane  away  from  the  focal  plane  of  a  single 
lens  system  introduces  a  2-D  wavefront  phase  error  that  can  be  parameterized  by  the  Zemike 
polynomial  for  defocus. 


To  produce  the  wave  model  of  a  telescope,  the  aberrations  are  summed  together  in 
its  exit  pupil,  where  the  aberration  free  pupil  transmittance  function,  A(u,v),  has  the 
coordinates  u  and  v  in  the  pupil  plane.  Wavefront  error,  W(u,v),  caused  by  defocus  is 
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introduced  into  the  pupil  function  using  the  Zernike  polynomial  for  defocus,  <t>4(u,v) , 


which  is  [31] 


fa(u,  v)  =  3.464(u2  +  V2)  - 1 .732. 


(2.1) 


The  amount  of  focus  error  is  captured  by  scaling  </>4(u,v)  with  a  Zemike  coefficient  for 
defocus,  Z4,  such  that 

W(u,v)  =  Z4  -</>4(u,v).  (2.2) 


An  image  of  the  unsealed  <p4(u,  v)  is  shown  at  the  bottom  of  Figure  4. 


2.1.2  Zernike  Polynomials  and  the  Generalized  Pupil  Function 
Compressing  the  notation  from  two  dimensions  (2-D)  to  one  (1-D)  for  simplified 
presentation,  the  wavefront  error,  W  (ux),  in  an  optical  system  can  be  decomposed  into  N- 

number  of  Zemike  polynomials,  (/\{ux)- (f>N{ux),  represented  as 

W(u1)  =  Zi  ■  (j>x  (uj)  + ...  +  ZN  -^N{ux),  (2.3) 

where  ux  is  a  coordinate  in  the  pupil  plane  and  Zx  -  ZN  are  the  Zernike  coefficients  [31]. 
Images  of  the  first  eleven  lower  order  Zemike  polynomials,  (j)2  to  <j)xx,  except  piston,  <j)x, 

are  shown  in  Figure  5.  Each  of  these  lower  order  polynomials  can  be  associated  with 
aberrations  that  arise  from  imperfection  in  the  optical  design  or  optical  alignment  and  are 
used  in  the  SST’s  current  alignment  process  [17]. 
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Figure  5.  Images  of  Zemike  polynomials  numbers  2-1 1 


The  aberrated  wavefront  error  is  represented  in  the  generalized  pupil  function, 
T(ut),  as 

T(ul)  =  A(ul)exp[j-W(ul)'],  (2.4) 

where  A(ut)  is  the  pupil  transmittance  function  [31].  From  the  generalized  pupil 

function,  a  model  of  the  electric  field,  H(m)  in  the  image  plane  as  a  function  of  image 
plane  pixel  coordinates,  m,  is  computed  using  a  discrete  Fourier  transfonn, 

H(m )  =  y  P( ux )  exp  ( jlnmu , ) 

(2.5) 

and  the  corresponding  optical  PSF,  h(m)  is 

h(m)  =  |//(m)|2  =  .  '  (2.6) 

2.1.3  Intensity  Model  and  PSF  Estimation 

The  images  of  stars  that  are  recorded  by  a  telescope  can  be  modeled  using  its  PSF. 
The  intensity  model  for  the  star  irradiance  centered  on  the  optical  axis  is 

i(m)  =  0X  •  h(m)  +  B.  (2-7) 
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It  includes  additional  terms  to  account  for  the  background  light,  B,  and  the  total  number 
of  photons  emitted  from  the  star  per  integration  time,  0X  [32].  Having  an  accurate  PSF 

and  intensity  model  for  the  stars  is  critical  for  both  the  phase  retrieval  and  detection 
algorithms  that  are  explored  in  this  work. 

2.2  Phase  Retrieval  from  Stellar  Images 

Phase  retrieval  is  a  set  of  techniques  used  to  determine  the  phase  aberrations  of  an 
optical  system  using  point  source  image  intensity  data  (i.e.  system  impulse  response). 

The  SST  phase  aberrations  must  be  recovered  using  post  processing  because  the  focal 
surface  array  cannot  measure  the  phase  directly.  The  existing  phase  retrieval  methods 
identified  from  literature  include  curvature  sensing,  least  squares  fitting  and  the 
Gerchberg-Saxton  algorithm  [33,  34,  35].  The  predominate  limitation  of  curvature 
sensing  and  least  squares  fits  is  that  they  require  the  point  source  image  to  be  defocused 
to  generate  enough  phase  diversity  to  distinguish  the  contribution  of  each  of  the 
polynomial  to  the  total  wavefront  error  and  corresponding  intensity  in  the  image  plane 
[17,  35],  The  Gerchberg-Saxton  (GS)  algorithm  on  the  other  hand  is  limited  because  it 
produces  a  wrapped  2-D  phase  for  larger  aberrations.  Attempts  have  been  made  to 
produce  reliable  2-D  phase  unwrapping  algorithms,  but  their  perfonnance  is  ultimately 
limited  by  branch  cuts  and  noise  [36,  37]. 

All  the  aforementioned  phase  retrieval  techniques  require  short  exposure  images 
to  work  in  the  presence  of  an  atmosphere.  When  long  exposure  imagery  is  used  for  phase 
retrieval  these  methods  prove  unreliable  [18].  A  method  that  works  to  some  extent  with 
the  long  exposure  imagery  and  has  been  used  to  estimate  the  SST’s  optical  prescription  is 
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the  Donut  software.  However,  the  phase  retrieval  with  Donut  was  only  conducted  on 
nights  with  excellent  seeing,  which  minimized  the  effects  of  the  atmosphere  [17]. 

The  current  process  for  alignment  of  the  SST’s  optics  is  critically  reliant  on  the 
Donut  software.  The  Donut  software’s  dependence  on  curvature  sensing  for  its  initial 
estimates  of  the  lower  order  Zemike  aberration  coefficients  (Zi-Ze)  is  a  problematic  point 
because  the  SST  alignment  process  requires  Zernike  coefficients  up  to  at  least  Z\\.  A 
second  weakness  in  the  method  is  that  curvature  sensing  does  not  account  for  the 
atmospheric  blurring,  so  poor  seeing  conditions  will  cause  error  in  the  initial  Zernike 
estimates.  To  mitigate  this  known  limitation,  the  site  engineers  monitor  the  seeing 
conditions  and  only  align  the  optics  on  nights  with  good  atmospheric  conditions 
(r0  >  10  cm) .  The  third  weak  point  in  the  algorithm  is  that  the  telescope  must  be  out  of 

focus  to  work,  but  in  the  process  of  refocusing  the  telescope  the  aberrations  are  likely  to 
change. 

With  future  SST  sites  planned  in  places  with  less  ideal  atmospheric  seeing 
conditions  than  New  Mexico,  there  will  be  a  need  to  correctly  estimate  the  Zernike 
coefficients  in  poor  atmospheric  conditions.  For  instance,  sites  surveyed  in  Australia 
have  seeing  parameter  estimated  as  poor  as  an  r0  of  6  cm  on  average.  Therefore,  finding 
a  better  phase  retrieval  method  that  can  more  reliably  estimate  higher  order  Zernike 
coefficients,  account  for  long  exposure  atmospheric  effects,  and  work  in-focus  should 
help  to  better  align  the  SST  and  produce  a  more  accurate  PSF  model. 

In  this  dissertation,  two  new  phase  retrieval  approaches  that  can  improve  the 
perfonnance  of  the  SST  are  described.  Both  techniques  achieve  joint  estimation  of  the 
static  telescope  aberrations  without  requiring  defocusing  of  the  telescope.  This  is  more 
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useful  because  moving  the  SST’s  secondary  and  tertiary  mirrors  for  defocusing  is  not 
practical  on  a  regular  basis.  In  addition,  changes  to  the  phase  errors  over  time  are  not 
measurable  with  the  current  phase  retrieval  techniques  from  focused  star  images  so  the 
phase  errors  cannot  be  continuously  monitored  or  removed.  The  estimation  of  the 
telescope’s  PSF  model  from  a  focused  spot  is  used  for  the  MHT.  Also,  a  more  accurate 
PSF  could  improve  the  SST’s  detection  perfonnance. 

2.2.1  CRLB  for  Zernike  Coefficients 

In  order  to  explore  the  theoretical  perfonnance  of  jointly  estimating  Zernike 
coefficients  for  unbiased  estimators,  the  Cramer-Rao  lower  bound  (CRLB)  is  used  to 
provide  the  statistical  lower  limit  for  variance  of  estimated  Zernike  coefficients  [38]. 
Fienup  et  al.  derived  the  CRLB  for  estimates  of  Zernike  coefficients  to  characterize  the 
aberrations  in  the  Hubble  Space  Telescope  (HST)  using  phase  retrieval  [32].  The  CRLB 
derived  for  the  HST  Zernike  estimates  did  not  include  a  parameter  for  the  atmosphere 
because  the  HST  is  above  the  atmosphere.  A  derivation  for  the  CRLB  for  Zernike 
coefficient  estimates  that  includes  the  seeing  parameter  in  Chapter  III  provides  a  tool  for 
comparing  the  theoretical  perfonnance  of  the  long  exposure  atmospheric  phase  retrieval 
method  with  phase  retrieval  when  no  atmosphere  is  present  [19]. 

2.2.2  Least  Squares  and  Gerchberg-Saxton  Phase  Retrieval 

In  their  paper,  Krist  et  al.  describe  a  Levenberg-Marquardt  least-squares  (LS) 
method  to  estimate  the  low  order  Zernike  coefficients  from  imagery  data  for  the 
aberrations  in  the  HST  before  and  after  conection  [34,  35].  In  order  for  the  LS  technique 


17 


to  work,  it  required  them  to  decouple  the  aberrations  by  defocusing  the  telescope  to 
introduce  phase  diversity.  In  Krist’s  paper,  the  lower  order  aberrations  are  determined 
and  the  remaining  phase  errors  are  retrieved  using  Gerchberg-Saxton  (GS)  phase  retrieval 
[34].  The  decomposition  of  the  remaining  phase  into  Zemike  polynomials  is  possible 
because  the  phase  errors  were  small  enough  that  they  remained  unwrapped.  Using  the 
orthonormality  of  the  Zernike  polynomials,  they  detennine  the  magnitude  of  the  Zernike 
coefficients  from  the  GS  phase. 

The  process  for  GS  phase  retrieval  is  illustrated  in  Figure  6  [34].  The  algorithm 
begins  with  initialized  wavefront  error,  W(ux),  by  multiplying  reasonable  initial 

conditions  for  the  Zernike  coefficients,  (Z1,Z2,...,ZAr),  with  their  respective  Zernike 

polynomials,  represented  as 

W(ux)  =  Z,  (ux )  +  ...  +  ZN  •  c/>N  (ux ) ,  (2.8) 

where  ux  is  a  coordinate  in  the  pupil  plane  and  N  is  the  total  number  of  Zernike 

coefficients  to  be  estimated.  The  aberrations  are  represented  in  the  generalized  pupil 
function  as 

T(ux)  =  A{ux)Qxy[j -W  {ux)\  (2.9) 

where,  A(ux)  is  the  pupil  function.  From  the  generalized  pupil  function  a  model  of  the 

electric  field,  H(m,W),  in  the  detector  plane  as  a  function  of  CCD  pixel  coordinates,  m, 
and  pupil  plane  wavefront,  W,  is  computed  using  a  discrete  Fourier  transform, 

H(m,W )  =  Y.PQO exp ( j'27imul )  =  yjh(m)  exp[  j  ■  xu  (/«)]  (2.10) 
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where  m(m)  is  the  phase  of  the  electric  field  in  the  image  plane.  The  square  root  of  the 


measured  PSF  estimated  from  the  image  intensity  data,  yjh  (m),  replaces  the  square  root 

of  the  predicted  PSF,  sjh  ( m  ) ,  which  is  output  from  the  fast  Fourier  transfonn  (FFT). 
Then  the  new  field  is  inversely  Fourier  transformed  to  obtain  a  new  estimate  for  the 
phase  in  the  pupil  plane.  Next,  the  estimated  pupil  function,  A'(u i),  is  replaced  by  the 
known  pupil  and  the  process  is  repeated  for  a  fixed  number  of  iterations  (300  iterations 
for  the  GS  phase  retrieval  algorithms  used  in  this  dissertation)  to  estimate  the  phase, 

)  [3 1  ].  While  the  technique  provides  a  wrapped  phase  that  is  difficult  to  unwrap  for 
Zernike  decomposition,  it  also  provides  estimates  of  the  electric  field  in  the  focal  plane. 
The  estimated  electric  field,  H(m),  proves  useful  in  estimating  the  Zernike  coefficients 
with  the  short  exposure  atmospheric  technique  described  in  Chapter  V. 
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Figure  6.  Gerchberg-Saxton  phase  retrieval  block  diagram 


2.2.3  Curvature  Sensing 

Roddier  developed  the  curvature  sensing  method  to  determine  the  low  order 
aberrations  in  a  telescope  by  de focusing  the  image  [33].  The  presence  of  aberrations 
causes  the  intensity  of  the  edge  of  the  defocus  spot  to  curve  in  a  way  that  can  be  used  to 
determine  the  Zernike  coefficients.  Later,  a  single  plane  technique  for  curvature  sensing 
was  developed  by  Hickson  [39].  Single  plane  curvature  sensing  gives  the  initial 
estimates  used  by  the  Donut  algorithm  [17].  Those  estimates  are  then  locally  adjusted  to 
refine  the  estimate  using  the  Zemax®  optical  ray  tracing  software  and  tested  on  the 
telescope  to  see  if  the  aberration  estimates  become  better  or  worse  in  an  iterative 
alignment  process. 
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2.3  Atmospheric  Models 

Two  different  models  for  the  atmosphere  can  be  used  to  determine  the 
atmospheric  effects  on  the  total  phase  and  PSF  of  the  optical  system.  The  long  exposure 
atmospheric  model  applies  when  exposure  times  of  images  viewed  through  the  Earth’s 
atmosphere  are  much  greater  than  10  ms  [20],  Whereas,  the  short  exposure  atmospheric 
model  holds  for  exposure  times  that  are  much  less  than  10  ms.  The  two  different  models 
affect  the  total  phase  and  PSF  in  different  ways,  and  therefore  change  the  method  of 
phase  retrieval  techniques  required  to  determine  an  accurate  model  of  the  telescope’s 
phase  errors  and  PSF. 

2.3.1  Long  Exposure  Atmosphere 

Because  SST  uses  a  shutter  with  an  integration  time  greater  than  25  ms,  an 
accepted  model  for  that  atmosphere  is  a  long-exposure  atmospheric  transfer  function, 
which  is  defined  by  Goodman  as  [20] 

exp<  -3.44  ^  f  Ul 
l  ro 

In  Eq.  (2.11),  A  is  the  mean  wavelength,  /  is  the  telescope  focal  length,  u2  is  spatial 
frequency,  and  r0  is  the  atmospheric  seeing  parameter.  The  transfer  function  of  the 
optical  system,  J-f  t  ( u2 ) ,  can  be  detennined  from  by  taking  the  Fourier  transform,  J7, 
of  the  modeled  optical  PSF,  hopt(m),  as 

av(ui  )=f[h,jmA  <212) 


W  atm  (U2) 
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The  long  exposure  PSF,  hL(m),  is  then  computed  using  the  inverse  Fourier  transfonn, 

'f  _1,  of  the  combined  atmospheric  and  optics  transfer  functions.  The  effect  of  the  finite 
square  pixel  width,  a,  is  detennined  using  a  rectangle  function  with  the  following  transfer 
function 


^ PiXei  («2  )  =  IF  [reel (a  •  m)] . 

The  modeled  PSF  centered  on  a  pixel  is  then  computed  as 

hL  M  =  J"‘  [Mopt  (M2  )  •  ^ pixel  (M2  )  •  rtatm  (u2  ) 


(2.13) 


(2.14) 


By  including  the  effects  of  the  long  exposure  atmosphere  in  the  PSF,  the  point  source 
image  intensity  model,  iL(m),  becomes 


iL(m)  =  0X  ■  hL(m )  +  B. 


(2.15) 


2.3.2  Short  Exposure  Atmosphere 

For  ground  based  telescopes  with  frame  transfer  cameras,  such  as  potential  future 
variants  of  the  SST,  the  exposure  time  can  be  limited  to  a  period  where  the  short 
exposure  model  applies.  The  short  exposure  atmosphere  can  be  described  by  a  set  of 
Zernike  polynomials,  which  have  coefficients  that  are  zero  mean  Gaussian  random 
variables  [27,  28,  40],  Therefore,  the  wavefront  error  for  a  telescope  with  a  short 
exposure  atmosphere,  Wtotal  s,  can  be  represented  by  expanding  Eq.  (2.3)  to  include  the 

atmospheric  contribution  to  the  Zernike  coefficients,  Z?  atm  to  ZN  atm ,  along  with  the 
static  telescope  optical  aberrations  Zernike  coefficient,  Z-,  opt  to  ZN  opt.  The  total  short 
exposure  wavefront  error,  Ws,  from  both  the  static  optical  aberrations  and  the 
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atmospheric  aberrations  can  be  describe  using  the  linear  combination  of  the  scaled 
polynomials 

W$  =  {^2 _opl  +  Z2_atm  )  '  ^2  +  •"  +  N _opt  +  ^  1 M  atm  )  '  ft N •  (2-16) 

By  using  the  combined  atmosphere  and  telescope  wavefront  error  from  Eq.  (2.16)  and 
then  inserting  it  into  Eqs.  (2.4),  (2.5)  and  (2.6)  the  combined  PSF  of  the  short  exposure 
atmosphere  and  the  telescope,  hs(m),  and  is  found  as, 

2 

hs(m)  =  £A(ui)-eJWsM-eJ2*mu'  .  (2.17) 

“i 

To  include  the  pixel  effect  the  short  exposure  OTF,  J-fs  (u2),  is 

^s(u2)  =  T[hs(mj],  (2.18) 

and  the  pixilated  short  exposure  PSF,  hs  pjx  (/;?),  is 

K_,t  (m)  =  T'1  |X  (« , )  •  («,)].  (2.19) 

Then  the  short  exposure  intensity  model,  is(m),  is  computed  as 

is(fn)  =  6>  •  hS  pix(m )  +  B.  (2.20) 

Accurate  modeling  of  the  long  and  short  exposure  PSF  and  corresponding  intensity 
models  is  critical  for  both  long  and  short  exposure  phase  retrieval.  In  addition,  a  MHT  is 
not  possible  without  an  accurate  long  exposure  image  intensity  model,  iL(m). 

2.4  Detectors 

Two  types  of  optical  detection  processes  of  unknown  space  objects  with  ground 
based  telescopes  are  discussed  in  the  literature.  The  first  type  of  detection  algorithm  is 
based  on  images  with  a  limited  exposure  time  such  that  objects  in  the  field  of  view  can  be 
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treated  as  point  sources  [25,  41,  42],  In  the  other  type  of  detection  algorithm,  the 
integration  time  is  longer  and  the  system  is  not  tracking  the  object  causing  each  of  the 
moving  objects  to  streak  over  multiple  pixels  within  the  frame  [23,  24],  The  image 
processing  technique  for  the  detection  of  space  objects  used  by  the  SST  limits  the 
exposure  time  such  that  the  objects  can  be  treated  as  point  sources.  Therefore,  the 
detection  algorithms  covered  in  the  comparative  analysis  in  this  body  of  work  are 
designed  to  detect  point  sources  not  streaks. 

The  three  single  frame  detection  schemes  based  on  point  sources  images 
identified  to  date  are  background  suppression  normalization  &  binary  quantization,  linear 
threshold  correlation,  and  non-linear  threshold  correlation  [25,  41,  42].  All  three 
detectors  can  be  derived  from  a  binary  hypothesis  test  expressed  as  the  following 
likelihood  ratio  tests  (LRT)  for  unifonn  cost  and  equal  priors  [1] 

_  P(d{w,z)V{w,z)<=[\,Md\\Hl)> 
P(d(w,z)\/(w,z)e[l,Md]\H0)  < 

Ho 

where  d(w,z)  is  the  image  data,  w  and  z  are  integer  pixel  coordinates  and  Md  is  the 
number  of  pixels  in  one  dimension  of  a  chosen  square  window  in  detector  plane.  In  this 
case,  //,  is  the  hypothesis  that  an  object  is  present  in  the  pixel  of  interest  and  H0  is  the 
hypothesis  that  an  object  is  not  present  in  the  pixel  of  interest.  The  joint  conditional 
probability  of  the  data  given  hypothesis  Hj,i  e  {0,1}  is  true,  is 

P(d(w,z)V(w,z)e[\,Md]\Hi ). 
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2.4.1  Point  Detection 


The  detector  currently  used  by  the  SST  is  based  on  background  suppression 
nonnalization  and  binary  quantization,  which  is  a  process  of  detecting  an  object  in  a 
single  CCD  pixel.  It  will  be  also  referred  to  as  either  a  point  detector  or  the  baseline 
detector.  As  stated  at  the  beginning  of  the  chapter,  the  detection  algorithm  currently 
used  in  the  SST  is  adopted  from  the  algorithm  developed  for  LINEAR  [25]. 
Mathematically,  the  single  frame  point  detection  of  an  object  from  the  SST  imagery  data, 

d(cx,cy  j,  in  a  pixel  with  coordinates  (■ cx,cy )  is  performed  as 


SV%,c,) 


(d(c  ,c  )- B)  > 

- ; -  r 

<j  < 

Ho 


(2.22) 


where  B  is  the  local  background,  a  is  the  standard  deviation  of  the  noise,  and  y  is  the 
detection  threshold.  The  background,  B ,  can  be  computed  as  the  local  sample  median  of 
the  data, 

B  =  median  [j  (w,z)  V(w,z)  e  [l,Mrf]],  (2.23) 


and  the  local  standard  deviation  is 


<j 


M  d  Md 


V=1  z= 1 


B 2, 


(2.24) 


where  Md  is  the  number  of  pixels  in  one  dimension  of  a  chosen  window  in  the  detector 


plane  centered  on  the  pixel  of  interest,  (yCx,cy  j.  Pixels  with  a  SNR  greater  than  the 

detection  threshold  are  classified  as  containing  a  target  and  passed  on  for  further 
processing. 


25 


For  this  method  of  single  frame  point  detection,  the  SNR  is  degraded  significantly 
because  the  SST’s  PSF  is  much  larger  than  the  size  of  a  single  pixel,  therefore  the  SST 
data  is  binned  into  2  by  2  pixels.  While  this  method  is  relatively  effective,  its 
perfonnance  is  still  inhibited  by  two  key  physical  limitations.  The  first  issue  with  this 
detector  is  that  the  telescope  cannot  focus  the  light  from  a  star  into  a  binned  pixel  causing 
a  decrease  in  SNR  output  of  the  point  detector.  The  second  is  that  the  objects  are  not 
always  centered  on  a  single  binned  pixel  so  that  when  an  object  falls  in  the  corner  or  side 
of  a  pixel  the  SNR  is  greatly  reduced. 


2.4.2  Correlation  versus  Point  Detection 

The  correlator  is  designed  to  achieve  a  chosen  probability  of  false  alarm,  PFA, 
under  the  H0  case  and  the  image  noise  is  modeled  as  Gaussian,  which  matches  the  SST 

noise  distribution.  A  Poisson  distribution  for  noise  would  also  be  equally  valid,  however 
Pohlig’s  derivation  using  that  assumption  led  to  a  detector  that  was  dependent  on  target 
irradiance  [42].  To  remove  the  detector’s  dependence  on  target  irradiance  a  log 
approximation  is  made  assuming  that  the  target  irradiance  is  low.  Then  the  paper 
concedes  that  the  distribution  of  the  noise  for  these  dim  objects  would  not  be  Poisson,  but 
have  a  similar  distribution. 

By  choosing  to  use  a  Gaussian  noise  distribution  the  LRT  becomes 


Ag 


P[d\Hx] 

P[d\H0\ 


Md  Md 


w= 1  z= 1  ~42kcj  > 


M.  M, 


1 


1, 


TT-TT  1  \  < 

nriTrr-^  f  »• 

W=1  Z=1  V  CJ 


(2.25) 
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where  w  and  z  are  pixel  locations  in  the  window,  and  the  long  exposure  PSF  is  hL  (w,  z ) . 
The  value  Md  is  the  total  number  of  pixels  in  the  window,  B  is  the  background  photon 

count  in  the  image,  0  is  the  space  object’s  irradiance,  and  a  is  the  standard  deviation  of 
the  noise.  The  sufficient  statistic  for  the  LRT  is  designed  to  maintain  the  same  false 
alann  rate  as  the  baseline  detector,  which  is  determined  by  the  H0  case.  Taking  the 
natural  log,  Eq.  (2.25)  reduces  to  the  following  fonn 


Md  Md 


w= i  z= i  2cr 


> 


log(AG)  =  J]£— T  -2B-6hL(w,z)  +  2d(w,z)-dhL(w,z)-(6hL(w,z)y  0.  (2.26) 


J< 


Since  the  PSF  can  be  estimated  independently  from  auxiliary  processes,  Eq.  (2.26)  can  be 
rearranged  as 


h , 


Yy£j(d(w,z)-B)hL(w,z)>  ^^Y\{hL(w,z)) 

w=  1  z= 1  ^  ^  vi=!  z-1 


<  2tr 

H , 


(2.27) 


The  selection  of  6  will  be  chosen  to  achieve  the  desired  threshold.  To  convert  Eq.  (2.27) 
into  a  sufficient  statistic  in  terms  of  signal-to-noise  ratio  (SNR),  the  background 
suppressed  data  is  a  new  random  variable,  d2,  with  zero  mean  in  the  H0  case  [2], 


d2(w,z)  =  d(w,z)~  B.  (2.28) 

The  correlation  of  the  PSF  with  the  background  suppressed  data  then  becomes 

Md  Md 

X  Z  d2  ( w’  z)hL  (w  - , z  -  s )’  (2-29) 

w—\  z—1 


where  cx  and  cy  are  the  coordinates  of  the  pixel  being  tested.  The  resulting  quantity  also 
has  a  mean,  //2,  where 
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Ml  =  E 


Mj  M, 


_ w=  1  z= 1 
Md  Mj 


(2.30) 


E  E  E  [ di  (  w,  z)]^  [ ^ ( w  ~  cx, z  ~  cy  )]  =  °, 


W=1  z= 1 


and  a  variance,  a of 


cj;=E 


=  E 


(  Md  Md 


E  E  ‘ d2  ( ' W’Z)hL  (  W  -  Cx ’  Z  -  Cy  ) 

V  W=1  z= 1 


Md  Md 


j 

Mj  Mj 


E  E 1 ^2  ( ■ w,  Z  K  ( W  -  cx  ,  z  -  cy )  E  E  ^2  ( m > «  K  (rn-cx,n-  cy  ) 

w=l  z=l  m=l  n=l 

.  Md  Mj  Mj 

=  E  E  E  E E  [d2  ( z)] ^  [d2  ( m,n)JiL  (w  -  cx,  z  -  cy)hL  (in  -  cx ,  n  -  cy). 


(2.31) 


w=l  z=l  m=l  «=1 

Eq.  (2.3 1)  can  be  simplified  using  two  cases,  one  when  w  ^  m  and/or  z^n  and  the  other 
when  w  =  m  and  z  =  n.  Using  the  Kronecker  delta  function,  S(w-m,z-n), 


Mj  Mj  Mj  Mj 


ct2  =  EEEE^  (w,z)]E[d2  ( m,n )] 

w=l  z= 1  m—l  n—l 

xh(w - cx,z  - cy)hL(m  -cx,n  -cy)(\-8(w-m,z-n)'} 

Md  Mj 

+EE£[^2  2{w,z)fil(w-cx,z-cy)S(w-m,z-n) 

W=1  z= 1 

Mj  Mj 

=cx2EE/ii2(w,z). 


(2.32) 


=1  Z— 1 


Therefore,  the  standard  deviation  of  the  normalized  data  convolved  with  the  total  system 
PSF  is 


\Mj  Mj 

C72=ffJZZ/f(W>Z)'  (2-33) 

\  VV=1  Z=1 

The  sufficient  statistic  for  the  correlator  is  then  found  by  dividing  Eq.  (2.27)  by  Eq. 

(2.33).  Then  the  LRT  reduces  to  the  following  correlation  operation  normalized  in  terms 
of  the  correlator’s  signal-to-noise  ratio,  SNRcorr, 
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M,  M, 


=  r-  (2.34) 


aJX  XhL2(W’Z) 


H „ 


aJ££hL2(w’z) 


W= 1  Z=1 


w=l  z=l 


The  results  of  the  sufficient  statistic  are  then  compared  against  the  selected  SNR 
threshold,  y,  which  is  set  to  achieve  a  desired  PFA. 

In  the  case  of  the  baseline  point  detector,  the  PSF  is  one  pixel  represented  as  a 
delta  function,  S(w  -  cx,z  -  cy).  The  sufficient  statistic  for  the  baseline  detector  in  terms 

of  the  point  detector’s  signal-to-noise  ratio,  SNRBaseline,  is 


M,,  Mj 


^^(c/(w,z)-£)£(w-Cv,z-cv) 


SNRBaseline  = 


C d(cx,cy)-B )  > 


M.,  Mj 


< 


Y 


(2.35) 


W=1  Z=1 


and  can  be  compared  against  the  same  threshold  as  the  correlator.  The  fact  that  both 
detectors  are  expressed  in  terms  of  SNR  and  use  the  same  threshold  makes  the 
comparison  of  the  two  detectors  possible.  This  is  because  for  each  pixel  being  tested  the 
detectors  will  produce  a  SNR  value.  When  comparing  the  two  detectors,  the  detector  that 
produces  the  higher  SNR  value  will  have  the  high  perfonnance. 


2.4.3  Undersampling  and  Correlation  Detection 

While  correlation  detectors  have  better  single  frame  detection  perfonnance  than 
the  point  detectors  currently  used  in  the  SST,  adequate  PSF  spatial  sampling  can  affect 
the  probability  of  detection,  Pd,  for  the  conelator  [25,  43].  O’Dell  et  al.  ’s  paper  shows 
the  effects  of  undersampling  images  on  the  performance  of  correlation  detectors  for  the 
optical  detection  of  space  objects.  What  O’Dell  doesn’t  describe  is  how  to  improve 
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detection  perfonnance  other  than  by  increasing  physical  spatial  sample  rates.  In  other 
words,  if  the  image  can’t  be  sampled  at  its  Nyquist  frequency,  can  the  aliasing  effect  be 
mitigated  using  a  modeled  PSF? 

To  understand  the  Nyquist  sampling  in  O’Dell’s  work,  objects  are  treated  as 
incoherent  point  sources  of  light  and  are  propagated  through  an  atmosphere  causing  the 
light  to  spread.  For  the  simulations  presented,  the  cutoff  frequency, 


2  *  Y 
f  =  0 

J  C  r>  9 

A  *  Z. 


(2.36) 


is  limited  by  the  atmospheric  seeing  parameter,  r0,  at  a  given  focal  length,  z.  and 
wavelength,  X  [31].  The  difference  between  sampling  to  meet  the  Rayleigh  criteria  is 
that  a  pixel  angle,  AR,  is 


A, 


1.22/1 


ro 


(2.37) 


and  the  Nyquist  criteria  has  pixel  angle,  A  v ,  of 


A, 


X 


2  -rn 


(2.38) 


The  curves  in  Figure  7  (a)  and  (b)  show  the  key  finding  of  O’Dell’s  paper  -  the 
correlator  has  a  higher  Pd  when  the  PSF  is  properly  sampled  [43].  When  the  object 
intensity  is  centered  on  a  pixel  the  Rayleigh  sampled  correlator  perfonnance  is  not  as 
drastically  degraded  as  compared  to  the  Nyquist  case  in  Figure  7  (a).  However,  when  the 
PSF  is  in  the  comer  of  a  pixel,  the  Rayleigh  sampled  correlator  has  a  significant  reduction 
in  detection  performance  as  shown  in  Figure  7  (b). 
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Figure  7.  Comparison  of  aliased  and  unaliased  detector  performance  for  a  fixed  probability  of 
false  alarm  (a)  Nyquist  and  Rayleigh  sampled  correlator  detection  performance  with  PSF  centered 
on  a  pixel,  (b)  Nyquist  and  Rayleigh  correlator  detection  performance  with  PSF  centered  on  the 
comer  of  a  pixel  [43], 


2.5  Conclusions 

The  mission  of  synoptically  searching  deep  space  for  unknown  space  objects  has 
driven  the  requirements  for  wide  FOV  three  mirror  telescopes  such  as  the  SST  and  the 
Large  Synoptic  Survey  Telescope  (LSST)  [1,21].  The  tertiary  mirror  makes  these 
complex  telescopes  considerably  more  difficult  to  align  and  focus  than  traditional  two 
mirror  telescopes  driving  the  need  of  better  phase  retrieval  methods  [21,  44].  To 
compensate  for  the  unavoidable  blurring  of  images  from  the  atmosphere  and  telescope 
aberrations,  larger  pixels  (or  binned  pixels)  are  used  to  increase  SNR  with  the 
unavoidable  consequence  of  causing  undersampling  of  the  data  [18,  43].  Detection  can  be 
improved  by  binning  pixels  to  add  together  the  spread  signal  with  the  SST’s  current 
detection  algorithm.  However,  this  approach  is  not  ideal  because  it  includes  the  addition 
of  read  noise  from  each  pixel.  Also,  the  PSF  may  not  be  centered  on  a  pixel  so  that  the 
irradiance  will  be  detected  across  multiple  binned  pixels.  Therefore,  a  detection 
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algorithm  is  needed  that  accounts  for  the  PSF  that  spans  more  than  one  pixel  and  spatially 
filters  the  image  with  a  model  of  the  PSF  to  perfonn  detection.  Review  of  current 
literature,  summarized  in  the  remainder  of  this  section,  has  revealed  that  existing  phase 
retrieval  algorithms  and  detectors  have  limitations  that  ultimately  limit  the  optical 
detection  of  space  objects. 

Phase  retrieval  of  Zemike  coefficients  has  been  successfully  used  to  align  the  SST 
in  good  atmospheric  conditions  (r0  >  10  cm)  by  defocus ing  the  telescope,  but  techniques 

that  work  in  less  ideal  conditions  and  in-focus  are  needed  to  maximize  telescope 
performance  [17].  A  phase  retrieval  technique  that  works  in-focus  with  any  atmospheric 
conditions  will  provide  the  necessary  aberration  infonnation  from  the  SST’s  standard 
imagery  data  for  focus  and  alignment.  This  can  be  accomplished  without  going  through 
the  complex  procedure  of  moving  the  secondary  and  tertiary  mirrors  to  defocus  the 
telescope  [18].  This  would  enable  diagnostic  monitoring  of  the  telescope  in  order  to 
maintain  focus  and  alignment  with  standard  imagery  data.  However,  the  current 
documented  phase  retrieval  methods  of  curvature  sensing,  the  Gerchberg-Saxton 
algorithm  and  least  squares  fitting  are  limited  in  their  ability  to  accurately  estimate 
Zernike  coefficients  with  focused  data  [33,  34,  35,  45],  Thus,  the  new  phase  retrieval 
methods  (discussed  in  Chapters  III  and  V)  that  work  in- focus  with  poor  atmospheric 
conditions  are  desirable.  In  addition,  the  phase  retrieval  algorithm  described  in  Chapter 
III  has  already  proven  useful  in  estimating  PSFs  for  inclusion  in  a  new  detection  strategy 
for  the  SST  covered  in  Chapter  IV. 

The  development  of  a  phase  retrieval  algorithm  that  works  in-focus  would  be 
useful  not  only  for  the  focus  and  alignment  of  ground  based  telescopes,  but  also  for  space 
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based  telescopes  like  the  James  Webb  Space  Telescope  (JWST)  projected  to  launch  in 
2018  [46].  Currently,  NASA  plans  to  use  a  phase  retrieval  algorithm  called  the  hybrid 
diversity  algorithm  (HD A)  that  is  based  on  the  Gerchberg-Saxton  algorithm  with  another 
process  to  perfonn  phase-unwraping  [36,  45,  47].  The  HDA  requires  phase  diversity 
similar  to  the  LM  least  squares  method  used  for  characterizing  the  HST  [32].  The 
problem  of  generating  phase  diversity  with  defocus  was  overcome  in  the  HST  by  moving 
the  focal  plane  [32],  however  due  to  JWST  three  mirror  design  and  segmented  primary 
mirror,  defocusing  is  not  as  simple  [44].  To  overcome  that  challenge,  the  JWST  has 
additional  optics  on  two  separate  wheels  that  can  be  rotated  into  the  optical  path  to 
generate  defocus  [47,  48].  The  JWST  could  potentially  remove  the  requirements  for 
these  additional  optical  elements  and  reduce  overall  program  risk  by  using  the  phase 
retrieval  algorithm  discussed  in  Chapter  V. 

In  addition  to  sub-optimal  phase  retrieval  methods,  the  SST’s  current  detection 
performance  is  limited  due  to  the  design  of  its  baseline  detector  [25].  A  correlator, 
similar  to  the  one  developed  for  pan-STARRS,  could  improve  the  perfonnance  of  the 
SST  over  the  baseline  detector,  but  it  is  ultimately  limited  by  the  SST’s  undersampled 
data  [16,  41,  43].  Other  correlation  methods  for  the  detection  of  space  objects  have  been 
developed  as  discussed  previously  in  this  chapter.  Those  methods  work  with  objects  that 
move  across  multiple  pixels  during  a  single  exposure  causing  streaks,  whereas  deep  space 
objects  imaged  by  the  SST  are  effectively  point  sources  due  to  short  integration  times 
[23,  24].  Each  of  the  aforementioned  detectors  are  based  on  a  Gaussian  parametric  model 
for  the  noise;  however,  one  other  detector  discussed  in  literature  was  based  on  a  Poisson 
distribution  of  the  noise  [33].  The  MHT  is  derived  from  a  Gaussian  LRT  in  Chapter  IV 
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and  works  with  point  source  data  to  improve  the  detection  performance  of  the  SST  over 
both  the  correlator  and  its  baseline  detector  by  compensating  for  both  blurring  and 
undersampled  data.  The  following  three  chapters  cover  the  research  accomplished  to 
improve  both  phase  retrieval  and  detection  performance  of  the  SST,  but  should  also  be 
extensible  to  other  astronomical  telescopes. 
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III.  Phase  Retrieval  with  a  Long  Exposure  Atmosphere 


One  key  image  processing  technique  that  improves  the  SST’s  detection  is 
modeling  the  telescope’s  wavefront  error  with  Zernike  polynomials.  Estimates  of  the 
Zernike  polynomial  coefficients  produced  by  phase  retrieval  are  used  to  better  focus  and 
align  the  telescope  immediately  improving  SNR  and  thus  the  detection  perfonnance  [17]. 
In  addition,  the  same  coefficients  are  used  to  model  the  telescope’s  point  spread  function 
(PSF),  which  can  be  used  to  improve  the  SST’s  detection  sensitivity  using  the  multi¬ 
hypothesis  test  discussed  in  Chapter  IV. 

The  critical  technology  that  enables  the  SST’s  6  deg  wide  field-of-view  (FOV) 
camera  (shown  in  Figure  8)  is  the  unique  set  of  curved  charged  coupled  device  (CCD) 
arrays.  The  set  of  12  separate  curved  CCD  arrays  are  tiled  together  in  a  6  by  2  mosaic 
that  fonn  a  surface  with  a  5  m  radius  of  curvature  that  alleviates  some  of  the  aberrations 
inherent  to  the  optical  design.  During  the  design  process,  a  choice  was  made  to  maintain 
the  wide  FOV  of  the  telescope  with  the  larger  format  (12288  by  8192)  mosaic  detector 
and  a  mechanical  shutter  rather  than  a  smaller  format  mosaic  detector  in  a  frame  transfer 
camera.  Since  the  current  camera  does  not  have  a  frame  transfer  capability,  a  high  speed 
shutter  was  developed  for  the  camera  with  a  minimum  exposure  time,  r,  of  25  ms.  This  is 
still  considerably  longer  than  the  r  <  10  ms  typically  associated  with  a  short  exposure 
image  [20]. 
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Sensor  Subsystem  (MIT/LL) 


CCD  and  Wide  Field  of  View  Camera  High  Speed  Shutter 


Figure  8.  The  SST’s  6  deg  wide  field-of-view  camera  and  high  speed  mechanical  shutter 

One  of  the  main  challenges  in  optimizing  the  SST’s  performance  is  to  reduce  the 
point  spread  function  (PSF)  through  focus  and  alignment.  The  pixels  in  the  CCD  are  15 
pm.  However,  the  pixels  can  be  two  by  two  binned  to  mitigate  the  degradation  in 
detection  sensitivity  of  the  point  detector  that  occurs  due  to  atmospheric  blurring  and 
telescope  aberrations.  Figure  9  shows  the  variations  of  full  width  half  maximum 
(FWHM)  blur  spot  in  terms  of  30  pm  binned  pixels  over  three  months  leading  up  to  the 
final  alignment  of  the  telescope.  The  FWHM  is  a  measurement  of  the  width  across  the 
irradiance  pattern  produced  by  the  telescope’s  image  of  a  point  source  at  half  the 
maximum  value. 
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Figure  9.  Variations  in  the  Full  Width  Flalf  Maximum  (FWHM)  of  the  SST’s  Point  Spread 
Function  (PSF)  measured  by  2  by  2  binned  15  pm  pixels,  (provided  by  MIT/LL) 


Ideally,  the  PSF  FWHM  would  be  within  one  30  pm  binned  pixel.  Reducing  the 
PSF  size  is  possible  by  accurately  determining  the  amount  of  focus  error  and  other 
aberrations  in  the  image  of  a  calibration  star,  then  adjusting  the  focus  and  alignment  to 
reduce  the  blur  spot  size.  However,  to  find  unbiased  estimates  of  the  telescope 
aberrations  from  the  star  image,  the  atmospheric  effects  must  be  accounted  for  in  the 
phase  retrieval  technique.  As  mentioned  in  Chapter  II,  the  engineers  only  aligned  the 
optics  on  nights  with  seeing  greater  than  10  cm.  With  future  SST  sites  being  surveyed  in 
Australia,  where  there  is  less  ideal  atmospheric  seeing  conditions  than  in  New  Mexico, 
there  will  be  a  need  to  correctly  estimate  the  Zernike  coefficients  in  poor  seeing 
conditions.  This  is  one  of  the  reasons  that  there  is  a  need  for  improved  phase  retrieval 
techniques. 
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3.1  Telescope  Model 

The  telescope  model  introduced  in  the  background  chapter  and  used  for  this 
analysis  follows  a  similar  model  employed  for  phase  retrieval  of  the  HST  aberrations 
before  and  after  correction  [32].  The  SST  is  considered  a  linear  shift  invariant  system 
with  the  impulse  response  of  the  system  being  the  PSF.  The  light  propagating  from  the 
distance  point  source,  d{x),  is  assumed  to  be  temporally  incoherent  in  the  image  plane 
with  coordinates,  x.  The  parameters  used  in  the  telescope  model  are  listed  in  Table  1 . 

Table  1.  Telescope  Model  Parameters 


Parameter 

Value 

Center  Wavelength 

500  nm 

Telescope  Pupil/Obscuration  Diameter 

3.5  m  /  1.8  m 

Telescope  Effective  Focal  Length 

3.5  m 

CCD  Pixel  Pitch 

15  pm 

Star  Irradiance  per  Frame 

~104  photons 

Background  Irradiance  per  Frame 

300  photons 

The  SST’s  pupil  transmittance  function,  A(u,  v),  is  defined  by  its  annular  aperture 

shown  in  Figure  10  (a),  where  U  and  V  are  coordinates  in  the  pupil  plane.  Wavefront 
error  caused  by  defocus  is  introduced  into  the  pupil  function  using  the  Zemike 
polynomial  for  defocus,  <f)4(u,v),  [31], 

<f>4(u,v)  =  3.464  •  ( u 2  +  v2)  - 1.732.  (3.1) 
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The  amount  of  focus  error  is  captured  by  scaling  </>4(u,v )  with  a  Zemike  coefficient  for 


defocus,  .Z4 .  This  product  is  used  to  express  the  wavefront  error,  W(u,  v),  as 


W(u,v)  =  Z4  -<f>4(u,v). 


(3.2) 


An  image  of  (f>4(u,v)  scaled  by  a  25  wave  coefficient  is  shown  in  Figure  10  (b). 
Compressing  the  notation  from  two  dimensions  (2-D)  to  one  (1-D)  for  simplified 
presentation,  the  aberrations  are  then  represented  with  the  pupil  plane  coordinates,  u] ,  by 
the  generalized  pupil  function,  T(ux),  as 


T(ux)  =  A(ux)Qxy[j  -W  (ux)\. 
Then  the  telescope’s  PSF,  hopt(m),  is  computed  as  [3] 


KPSm)  = 


(3.3) 


(3.4) 


where  m  is  a  pixel  coordinate  in  the  detector  plane.  An  image  of  hopt(m)  with  25  waves  of 
defocus  is  shown  in  Figure  10  (c).  The  large  amount  of  defocus  causes  the  PSF  to  have 


an  annular  shape. 


Figure  10.  Telescope  Model  (a)  Pupil  function  used  to  model  the  SST  (b)  Zemike  polynomial  for 
defocus  with  Z4  =  25  waves  (c)  Telescopes  PSF  with  Z4  =  25  waves  of  focus  error. 
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For  a  more  complete  model  of  the  telescope,  the  effect  of  the  finite  square  pixels, 
o=15  pm,  is  included  in  the  PSF  where  the  transfer  function  for  the  pixels,  3~C pixel  (u2), 

and  telescope,  J~Copt  ( u2 ) ,  respectively  are  represented  as  the  following  discrete  Fourier 

transforms,  J~ , 

d~CplxeI{u 2)  =  T[rect(a-m)],  and  (3.5) 

3-C„p,  (m2)  =  T[_hopt(m)],  (3.6) 


where  u2  is  the  spatial  frequency.  Then  the  PSF  for  the  telescope,  hlelescope(m),  can  be 
computed  via  the  transfer  functions  as 


Klescopei ™)  =  T'  \Mopt  (M2  )  '  ^ puel  {U2  ) 


(3.7) 


Bright  star  images  observed  by  the  SST  have  been  measured  to  be  shot  noise  dominated 
(see  Appendix),  so  the  image  data,  d (m),  is  considered  to  be  Poisson  and  has  a  mean 
value  that  is  equal  to  the  irradiance  of  light  in  that  pixel  [32] 

is[d(m)]  =  z(m).  (3.8) 

The  model  for  the  star  irradiance  centered  on  the  optical  axis  is 


1  telescoped)  =  8 W  '  Klescopei™  ~  X)  +  B 

X  W  • 

=  dl-K,escoPei™)  +  B- 

It  includes  additional  terms  to  account  for  the  background  light,  B,  and  the  total  photons 
emitted  from  the  star  per  integration  time,  6X .  Assuming  statistical  independence 

between  pixels,  the  joint  distribution  of  the  image  data  is  represented  by  the  Poisson 
probability  mass  function  and  details  of  the  choice  of  this  PMF  can  be  found  in  the 
Appendix, 
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iiwhi 


e  ilelesco^m)  •  i 


telescope 


(m) 


d(n 


d  ( m ) ! 


The  associated  log  likelihood,  Z(Z4),  equation  is 


f  \ 

L{Z4)  =  In  Y\P[d(m)] 

V  m  ) 


=  E{-W-)  +  <*(w)  ‘  lnWo/Je(^)  -  In  J(m)l}. 


(3.10) 


(3-11) 


3.2  Cramer-Rao  Lower  Bounds  (CRLB)  for  Variance 

The  CRLB  for  variance  presented  in  this  section  was  derived  in  order  to  evaluate 
the  performance  of  the  joint  estimator  (i.e.  phase  retrieval)  for  the  Zemike  coefficients 
and  the  atmospheric  seeing  parameter.  The  CRLB  for  estimates  of  Zernike  coefficients 
was  previously  derived  for  evaluating  the  phase  retrieval  performance  on  the  HST  [32], 
The  difference  between  the  SST  and  the  HST  is  that  the  HST  does  not  have  to  image 
through  the  earth’s  atmosphere.  The  CRLB  herein  provides  a  theoretical  lower  limit  of 
variance  for  unbiased  estimates  of  the  Zemike  coefficient  for  defocus,  Z4,  and  the 
atmospheric  seeing  parameter,  r0.  Simulations  presented  in  this  chapter  demonstrate  that 

phase  retrieval  using  a  least  squares  estimator  produces  unbiased  estimates  of  those  two 
parameters  [32]. 

The  bounds  in  Figure  12  illustrate  that  standard  deviations  of  Z4  on  the  order  of 

10'1  waves  are  possible  at  practical  light  levels  even  in  the  presence  of  long  exposure 
atmosphere.  To  detennine  the  CRLB  for  estimates  of  the  Zemike  coefficient  for  defocus 
the  Fisher  information,  J{ZA),  is  computed  via  the  following  calculation  [38], 
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J(Z4)  =  -E 


d2L(Z4 ) 


dZA 


(3.12) 


In  this  equation  the  CRLB  for  the  variance  of  Z4  is  defined  as 

var(Z4)  >  J  (Z4)~‘ .  (3.13) 

The  first  and  second  derivative  of  the  log  likelihood  function,  Eq.  (3.11),  respectively  are 

d  ( m )  •  6X 


8L(Z4)  __  y, 

dZA 


i telescope  ( 


9 


dkte, escape  W_,md 


8Z, 


d2L(Z4 ) 
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d(m)Ox  e 

d\lescope{m) 

d  ( m )  •  6X 

dK,escope{m) 

i telescope  (  ^0 

dZ24 

1 telescope  (W) 

1 

N 

co 

_ i 

(3.14) 


.  (3.15) 


The  resulting  Fisher  infonnation  is 
J(Zt)  =  -E 


's2l(z ,y 

-  V 

dh,elescope(m ) 

dz  2 

_  Zj 

m 

1 

& 

o 

So 

•Si 

_ 1 

az4 

(3.16) 


The  derivative  of  the  PSF  with  respect  defocus  is 


dh  telescoped )  _  rr-l 


dZA 


dH  (m)  .  8H  (m)  .  * 

//(/»)+  '  -H(m) 


dZ4  v  7  sz4 


(3-17) 


The  derivative  of  the  wavefront  in  the  detector  plane  with  respect  to  (w.r.t)  Z4  is 

dH{m ) 


8Za 


■  =  j  •  (u,)A (u,)eJzMui)eJ'2mnUl 


(3.18) 


=  j-  T{^{ux)A{ui) 


,7Z4A(«i) 


Thus,  recalling  that  for  arbitrary  variables  a  and  b\ 

j[(a  +  jb)  -  (a  -  jb )]  =  -2  •  b  =  -2  ■  Im(a  +  y'b) 
leading  to  the  first  derivative  of  the  PSF  to  be 
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where  «3  represents  the  pupil  coordinates.  The  resulting  Fisher  infonnation  for  the 
optical  system  containing  focus  error  is 

Az,)  =  4ll— l(^Tj'(Imfj[A(«J)^(«))-^*l-)]-W(m)'})-jr^,(U!)ir.(3.20) 

mVz*)y  L  v  /  j; 


As  discussed  in  the  background  chapter,  because  the  SST  uses  a  shutter  with  an 
integration  time  greater  than  25  ms,  an  accepted  model  for  that  atmosphere  is  a  long- 
exposure  atmospheric  transfer  function,  which  is  given  by  Goodman  as  [20], 


In  Eq.  (3.21)  A  is  the  mean  wavelength,  /  is  the  telescope  focal  length,  and  r0  is  the 
atmospheric  seeing  parameter.  The  long  exposure  PSF  is  then  computed  as, 

hlong(m)  =  J" 1  [Mopt  W  •  M pixel  (Ul)-rtatm  (w2)].  (3.22) 

Samples  of  the  three  transfer  functions  are  shown  in  Figure  1 1  (a-c)  to  illustrate 
how  the  pixels  and  atmosphere  reduce  spatial  frequency  content  of  the  diffraction  limited 
telescope’s  optical  transfer  function.  The  horizontal  axis  is  shown  in  terms  of  the  spatial 
frequency,  u2,  divided  by  the  cutoff  frequency,  u0,  for  the  annular  telescope  pupil 
function.  Reducing  the  spatial  frequency  of  the  optical  system  causes  the  PSF  to  broaden 
due  to  the  Fourier  transform  relationship.  As  the  focus  error  increases,  J-f  (u2  ) 

begins  to  increasingly  limit  the  spatial  resolution  of  the  telescope. 
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Figure  11.  System  model  transfer  function  examples,  (a)  Telescope  models  optical  transfer 
function  (OTF)  with  Z4  =  0.  (b)  15  pm  pixels  transfer  function,  (c)  Atmospheric  transfer 
function  with  r0  =  8  cm. 


By  including  the  effects  of  the  atmosphere  in  the  PSF,  the  image  intensity  model  in  Eq. 
(3.9)  becomes 


hong(m)  =  S(xXng(m  ~x)  +  B. 

x 

The  elements  of  the  Fisher  information  matrix,  /, 


J(Z4)  J{Z4,r0j 
J(Z4,r0)  J(r0) 


(3.23) 


(3.24) 


are  calculated  in  order  to  determine  the  CRLB  for  variance  of  Z4  in  the  presence  of  an 
average  atmosphere.  Using  the  log  likelihood  function  in  Eq.  (3.11)  and  taking  the 
second  derivative  of  Eq.  (3.14)  w.r.t.  r0  &  Z4, 
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,(3.25) 


and  (3.26) 
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(3.27) 
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Because  E^d  (w)J  =  ilong  (m)  the  elements  of  the  Fisher  infonnation  matrix  are 

J(Z*,r0)  =  —E 
j(r0)  =  ~E 
j{z4)  =  -e 

Then  the  derivatives  of  the  PSF  are 

dhio„g(m ) 
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(3.28) 
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(3.30) 
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,  and  (3.31) 
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The  CRLB  is  computed  by  inverting  the  Fisher  infonnation  matrix 

>r\ 


var  (Z4)  cov(  Z4,r0) 

cov(Z4,r0)  var(r0) 


(3.32) 


The  resulting  bound  of  the  standard  deviation  for  estimates  of  the  Zernike 
coefficient  for  defocus,  aLB  (Z4,r0)  are  plotted  in  Figure  12  for  cases  with  and  without  an 


average  atmosphere  present.  Atmospheric  turbulence  increases  the  bound  and  as  r0 
decreases,  the  effect  of  the  atmosphere  on  the  bound  increases  because  a  decreasing  r0 
represents  a  more  turbulent  atmosphere.  In  addition,  as  Z4  decreases  the  lower  bound 
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increases.  Therefore,  estimation  of  Z4  should  become  more  inaccurate  as  the  amount  of 


defocus  decreases. 


Figure  12.  A  plot  of  the  CRLB  for  the  standard  deviation  of  the  Zernike  coefficient  for  defocus, 
crLB  (Z4,r0),  for  a  range  of  Z4  values.  The  standard  deviations  are  shown  for  cases  with  no 
atmosphere  in  the  model  and  increasing  atmospheric  seeing  by  changing  r0. 

3.3  Parameter  Estimation 

The  method  of  least  squares  (LS)  estimation  was  used  to  estimate  the  Zernike 
coefficient  for  defocus,  Z4,  from  simulated  star  data.  The  LS  method  is  used  because  of 
computer  precision  challenges  encountered  in  the  maximum  likelihood  estimation 
approach,  particularly  in  evaluating  the  natural  logarithm.  The  primary  contribution 
comes  from  to  the  large  background  levels  inside  the  log-likelihood  function.  To 
accurately  estimate  the  defocus  parameter  from  the  simulated  star  data,  an  estimate  of  the 
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object’s  irradiance,  6X,  must  be  calculated  from  the  data  by  taking  the  derivative  of  Eqs. 
(3.34)  and  (3.35)  then  setting  them  equal  to  zero  to  get  the  generalized  function 


'y\d(m)  - B\- h(m) 

m _ 

2>2M 


(3.33) 


The  intensity  models  from  Eqs.  (3.9)  and  (3.23)  are  used  to  define  the  elements  of 
the  sum  of  squares  matrices  for  the  telescope  model  without  an  atmosphere, 

Qteiescope  (Z4),  and  with  an  atmosphere,  Q  (Z4,r0).  The  elements  of  each  matrix  are 
calculated  respectively  as 


Qteiescope  (^4  ) 


d telescope  (tn ’  ^ 4)  I  dlS(x)htelesCOpe(x-m)-B 


and 


Q/ong(Z 4,?o)  ^ 


dlong(m’Z 4^o) -Y.^S(XdhlonJX  ~  "0  “  B 


(3.34) 

(3.35) 


where,  dtelescope(m,Z4)  and  dlong(m,Z4,r0)  are  data  from  a  simulated  star  and  correspond  to 
a  grid  formed  of  possible  Z4  &  r0  values.  For  the  single  parameter  estimate  of  Z4,  a 
vector  of  values  of  Qtdescope  are  formed  with  each  element  in  the  vector  corresponding  to 

the  following  range  of  defocus  parameters,  Z4  =  0,  .25,  ...,29.75, 30  waves.  For  the  joint 
estimation  Z4  and  r0  ,  a  matrix  of  values  of  Q,  are  computed  with  input  parameters 

from  the  sets  Z4  =  0,  .25, ...,29.75, 30  waves  and  r0  =  2, 2.1, ...,9.9,10  cm.  Then  Z4  is 
determined  without  accounting  for  the  atmosphere  by  the  single  parameter  estimate, 

Z4  =  arg  min  (Qtelescope) ,  (3.36) 
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or  by  accounting  for  the  atmosphere  with  the  joint  estimator 


=  arg min  Q,  .  (3.37) 

Z4,r0  \=5:/ 

By  finding  Z4  for  multiple  image  frames  of  simulated  star  data,  the  results  are  then  used 
to  determine  the  sample  mean  and  variance  for  the  LS  estimator. 

3.4  Phase  Retrieval  Simulations 

Simulated  star  data,  d,  (m),  viewed  through  a  long  exposure  atmosphere  and 

defocused  telescope  are  modeled  in  order  to  evaluate  the  performance  of  the  LS  phase 
retrieval  technique.  Stars  are  simulated  as  system  impulses,  §{x),  and  then  the  effects  of 
the  atmosphere,  telescope,  defocus,  pixilation,  background  light  and  star  intensity  are 
introduced  using  Eq.  (3.23).  The  Nyquist  sampled  images  of  a  star  with  and  without 
atmospheric  effects  are  shown  in  Figure  13.  (a)  &  (b)  respectively.  The  pixilated  images 
of  those  same  stars  are  picture  in  Figure  13.  (c)  &  (d).  Shot  noise  is  simulated  in  the  star 
data  using  a  Poisson  random  number  generator. 
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Figure  13.  Simulated  stars  (a)  simulated  analog  image  of  a  star  with  Z4  =  18  waves  and  no 
atmosphere,  (b)  Simulated  analog  image  of  a  star  with  Z4  =  18  waves  and  an  average  atmosphere 
where  r0  =  8  cm  (c)  Simulated  digital  image  of  a  star  with  Z4  =  18  waves  and  no  atmosphere,  (d) 
Simulated  digital  image  of  a  star  with  Z4  =  18  waves  and  an  average  atmosphere  where  r0  =  8  cm. 


To  produce  the  plot  in  Figure  14,  dlong  (m)  is  generated  without  shot  noise  and 

with  focus  errors  ranging  from  3-24  waves  in  order  to  evaluate  the  LS  estimators  biases. 
Estimates  of  defocus  using  Eq.  (3.36)  are  made  on  simulated  stars  with  and  without  an 
average  atmosphere  present.  The  graph  shows  that  when  the  simulated  star  data  does  not 
have  an  average  atmosphere,  the  single  parameter  estimator  is  unbiased.  Also,  when  the 
average  atmosphere  is  introduced  to  the  data  the  single  parameter  estimator  has  a  defocus 
dependent  bias.  In  contrast,  the  results  of  the  joint  estimator,  Eq.  (3.37),  detennined  from 
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the  same  simulated  star  data  with  an  average  atmosphere  present  does  not  have  a 
significant  bias. 
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Figure  14.  Estimated  defocus  parameter  determined  from  simulated  star  data  with  no  noise 
present  to  investigate  phase  retrieval  biases.  The  blue  X  marks  the  single  parameter  estimate  for 
defocus  without  an  atmosphere  present  in  the  simulated  star  data.  The  green  circles  are  the  single 
parameter  estimates  of  Z4  where  r0  =  8  cm  in  the  simulated  star  data.  The  red  boxes  are  the  joint 
parameter  estimate  of  Z4  where  r0  =  8  cm  in  the  simulated  star  data. 


To  further  evaluate  the  phase  retrieval  performance  of  the  LS  joint  estimator,  shot 
noise  is  added  d,  (m)  to  form  multiple  frames  of  simulated  star  data.  Then  the 

estimator’s  results  mean,  is[z4],  and  standard  deviation,  as,  for  each  defocus  value  are 

plotted  in  Figure  15.  As  the  amount  of  defocus  in  the  simulated  stars  decreases  the 
standard  deviation  of  the  estimates  increases  significantly  due  to  the  narrowing  of  the 
PSF,  measured  as  the  FWHM  on  the  right  hand  side  of  the  plot.  As  the  PSF  narrows,  less 
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of  the  shape  of  the  PSF  can  be  discerned  from  the  star  images  increasing  the  standard 
deviation  of  estimates  of  the  defocus  parameter. 
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Figure  15.  The  joint  parameter  estimates  from  the  simulated  stars  with  shot  noise. 


and<T 


are  represented  as  blue  dots  and  error  bars  then  plotted  as  a  function  of  the  simulated  defocus. 
The  FWHM  of  the  PSFs  are  plotted  with  the  green  asterisks  as  a  function  of  the  same  simulated 
defocus. 


In  Figure  16,  the  same  sample  standard  deviation,  <js,  from  the  joint  parameter 
estimates  in  Figure  15  are  plotted  along  with  their  associated  CRLB,  aLB  (Z4,r0).  The 
aLB  ( Z4,r0 )  is  not  achieved  by  crv ,  but  the  standard  deviation  is  below  a  wavelength  until 

the  blur  spot  becomes  too  small.  Overall  the  joint  estimator  performs  well  for  phase 
retrieval  of  the  defocus  and  seeing  parameters  because  the  estimator  is  unbiased  and  its 
variance  is  manageable  since  many  frames  of  data  can  be  recorded  to  reduce  the  variance 
to  the  desired  range. 
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Figure  16.  Both  the  sample  standard  deviation,  <JS,  of  the  joint  parameter  estimates  (blue  dots) 
and  the  CRLB,  crLB  (Z4,r0),  (red  stars)  are  plotted  as  a  function  of  simulated  star  defocus. 

3.5  Laboratory  Demonstration 

A  demonstration,  illustrated  in  Figure  17,  was  conducted  to  show  that  the  LS 
estimator  for  Z4  and  r0  works  beyond  the  pristine  conditions  of  simulation.  The  setup 

includes  a  point  source  created  with  a  red  LED  and  a  pin  hole  to  create  a  point  source. 

The  light  emanating  from  the  pin  hole  can  be  considered  a  spherical  wave.  The  light  then 
propagates  to  the  2  mm  aperture  of  an  intentionally  defocused  camera.  From  the  single 
lens  forming  the  aperture,  the  light  is  imaged  onto  a  CCD  array  with  a  pixel  pitch  of  16 
pm.  The  data  is  recorded  and  the  joint  estimate  of  Z4  and  ro  is  made  using  the  LS 
estimator.  The  same  setup  is  used  for  the  second  half  of  the  demonstration,  only  this  time 
a  thermal  source  is  placed  in  front  of  the  aperture  to  simulate  a  turbulent  atmosphere. 
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The  data  is  recorded  and  again  the  joint  estimate  of  Z\  and  ro  is  made  using  the  LS 
estimator. 

The  demonstration  results  are  shown  at  the  bottom  of  Figure  17.  On  the  lower 
left  is  the  image  of  the  point  source  without  the  thermal  source.  The  LS  estimate  of  that 
blur  spot  has  a  Z4  =  .7895  waves  and  an  r0  =  2.27  cm.  Then  the  blur  spot  on  the  lower 
right  is  recorded  with  the  thermal  source  on  and  the  estimates  changed  to  a 
Z4  =  .4737  waves  and  a  r0  =  .018  cm.  While  truth  data  for  the  defocus  and  the  seeing 
parameters  was  not  available  for  this  demonstration,  the  values  are  consist  with  the 
observed  blurring  of  the  point  source  image.  The  results  from  laboratory  data  are  an 
indication  that  the  LS  estimator  is  working.  With  the  heat  source  on,  r0  drops  below  the 

diameter  of  the  aperture  causing  a  measurable  blur  in  the  image  of  the  point  source  as 
anticipated. 


Point  Source 

-  -• - 


■> 


Plane  Wave 
yL=650nm 


2 mm  ccd 

|  Ax  =  16//m 

I - <1 . 

Thermal  f  _  ,  ()cm 
Source 


1  Frame  LS  Estimates 


RstSoiiEKAigAtai 


Z4  =  .4737  waves 


r0  =  1 .8  mm 


Figure  17.  Phase  retrieval  demonstration  setup  and  results 
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3.6  Joint  Estimation  of  Spherical  Error,  Defocus,  and  Atmospheric  Seeing 

Ray  tracing  analysis  of  the  SST  optical  design  predicts  it  to  have  spherical  error. 
Therefore,  it  is  important  to  investigate  the  effects  of  the  presence  of  spherical  error  when 
estimating  the  Zemike  coefficient  for  defocus.  Once  again  stars  were  simulated,  however 
this  time  the  Zemike  coefficient  for  spherical  error,  Z, , ,  was  included  in  the  simulated 

star  data.  Then  estimates  for  Zn  were  made  by  extending  the  joint  estimator  to 


=  arg  min 


(3.38) 


The  plot  in  Figure  18  shows  that  jointly  estimating  the  atmospheric  seeing 
parameter  plus  the  Zemike  coefficients  for  spherical  and  defocus  produces  an  unbiased 
estimate  when  those  aberrations  are  present  in  the  optical  system.  However,  the  omission 
of  either  spherical  or  defocus  parameters  from  the  estimate  produces  a  bias  estimate  for 
the  parameter  of  interest. 
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Figure  18.  Defocus  and  spherical  error  estimation  results.  Estimated  defocus  and  spherical  error 
parameters  determined  from  simulated  star  data  that  has  both  defocus  and  spherical  aberrations 

through  an  average  atmosphere  with  no  noise  present.  The  red  X  marks  the  joint  Z4,  Zn,  r0 
parameter  estimate  for  defocus.  The  green  asterisk  are  the  estimates  of  Z4  using  only  the  joint 
Z4,  r0  estimator.  The  light  blue  box  marks  the  joint  Z4,  Zn,  r0  parameter  estimate  for  spherical 
error.  The  blue  circles  are  the  estimates  of  Zu  using  only  the  joint  Z4,  r0  estimator.  The 
estimator  achieves  the  correct  value  when  all  three  parameters  are  jointly  estimated. 

3.7  Conclusions 

Based  on  the  simulations  discussed  in  this  chapter,  phase  retrieval  of  the 
atmospheric  seeing  parameter  and  Zernike  coefficient  for  defocus  can  be  accurately 
detennined  as  long  as  they  are  jointly  estimated.  The  long  exposure  atmosphere  can 
affect  estimation  of  Z4  and  Zn  thus  demonstrating  the  need  to  account  for  r0  when 
estimating  the  SST’s  aberrations.  That  being  said,  ideally  an  estimator  that  includes 
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more  Zernike  terms  would  be  preferable,  so  that  the  optical  system  can  be  better 
optimized  and  the  PSF  better  characterized.  To  that  end,  Chapter  V  on  short  exposure 
imagery  demonstrates  a  phase  retrieval  method  of  more  Zemike  coefficients  with  a  short 
exposure  atmosphere.  Unfortunately,  the  grid  search  method  is  too  computationally 
burdensome  for  the  estimation  of  more  coefficients  and  so  a  direct  search  method  is  used. 

In  the  next  chapter,  the  grid  search  method  presented  in  this  chapter  is  used  for 
phase  retrieval  of  the  defocus  and  seeing  parameters.  This  method  was  possible  because 
the  SST’s  other  aberrations  have  been  previously  estimated  using  the  Donut  software,  so 
that  they  could  be  included  in  the  PSF  model  to  avoid  biasing  of  the  seeing  and  defocus 
parameter  estimates.  The  SST’s  aberrations  were  phase  retrieved  with  Donut  within  days 
of  the  collection  of  experimental  data  so  that  the  PSF  was  relatively  current.  By 
including  the  other  parameters  in  the  PSF  model  and  phase  retrieving  the  atmospheric 
seeing  parameter  and  defocus  parameter  that  change  on  a  temporal  basis,  an  up-to-date 
model  of  the  PSF  is  formed.  The  PSF  model  retrieved  using  this  combined  parameter 
estimation  technique  is  used  in  the  MHT  and  has  an  advantage  over  using  a  star  in  the 
FOV  for  the  PSF  because  the  modeled  PSF  can  be  shifted  without  aliasing.  However, 
since  the  PSF  may  change  over  time  it  would  be  preferable  if  all  the  Zemike  coefficients 
used  in  the  PSF  model  could  be  characterized  from  focused  SST  data  as  demonstrated  in 
Chapter  V. 
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IV.  Improving  Detection  using  Multi-hypothesis  Testing 


The  mission  requirements  for  the  SST  differ  from  typical  astronomical  telescopes. 
The  SST  is  designed  to  scan  deep  space  to  detect  unknown  space  objects  and  correlate 
their  orbits  rather  than  dwell  on  stellar  objects  over  relatively  long  periods  of  time  in 
order  to  characterize  them  [3].  In  this  sense,  the  SST  is  a  precursor  to  other  wide  field  of 
view  synoptic  search  programs  like  the  LSST.  Two  trajectory  matched  filter  approaches 
for  asteroid  detection  are  described  in  [23]  and  [24];  however,  these  particular  approaches 
are  not  further  investigated  in  this  work.  This  chapter  discuses  match  filtering  the  spatial 
shape  of  the  object  in  a  single  observation  similar  to  a  correlator  or  the  multi-hypothesis 
test  (MHT). 

4.1  Introduction 

Currently  the  SST  uses  an  algorithm  developed  from  a  binary  hypothesis  test 
(BHT)  to  detect  space  objects  in  a  single  image  [25].  The  two  hypotheses  are  1)  the  null 
hypothesis  that  a  space  object’s  image  is  not  in  a  pixel  (Ho)  and  2)  the  alternative 
hypothesis  that  the  image  is  in  a  pixel  (H\).  In  contrast,  a  MHT  is  proposed  for  single 
frame  detection  that  selects  from  the  hypotheses  that  the  image  is  in  the  center,  a  corner, 
or  a  side  of  a  pixel  (Hi-Hg)  in  addition  to  Ho.  Although  the  use  of  more  hypotheses  might 
increase  the  detection  perfonnance  of  this  scheme,  a  finite  number  of  hypotheses  must  be 
chosen  in  order  to  make  the  use  of  the  test  numerically  tractable.  The  results 
demonstrated  using  nine  alternative  positional  hypotheses  serve  to  demonstrate  the  utility 
of  the  MHT  over  a  BHT,  but  do  not  necessarily  represent  the  optimal  performance 
achievable  by  a  MHT.  Since  the  SST  sensor  is  dominated  by  readout  noise  rather  than 
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shot  noise  at  low  light  levels  these  tests  consider  the  noise  to  be  Gaussian;  however, 
another  method  of  BHT  developed  by  Pohlig  has  been  derived  using  a  Poisson  noise 
distribution,  which  is  not  used  for  comparison  in  this  study  [42]. 

The  experiment  conducted  as  part  of  this  research  is  designed  to  detennine  which 
type  of  algorithm  is  best  at  detecting  dim,  unresolvable  objects  in  space  on  a  single  frame 
basis.  Since  single  frame  detection  decisions  are  typically  used  as  input  to  multi-frame 
detection  and  tracking  algorithms,  a  superior  single  frame  detector  will  enhance  the 
perfonnance  of  any  synoptic  search  telescope  looking  for  NEAs  or  space  debris  using 
this  three  frame  coincidence  approach  [25].  In  order  to  perform  this  study,  we  chose  to 
observe  a  satellite  in  GEO  that  is  gradually  going  into  eclipse  behind  the  earth.  In  this 
scenario,  the  unresolvable  satellite  body  experiences  an  ever  decreasing  amount  of  solar 
illumination,  providing  a  range  of  intensity  values  over  which  to  test  the  performance  of 
different  algorithms.  Since  the  presence  and  location  of  the  satellite  is  simple  to  establish 
when  it  is  brightly  lit,  all  detection  algorithms  will  successfully  detect  the  object  before  it 
begins  to  go  into  eclipse.  The  telescope  is  pointed  directly  at  the  satellite  and  then 
observes  it  as  it  goes  into  the  shadow  of  the  earth.  Because  the  presence  of  the  object  is 
known  (and  further  verified  when  it  emerges  from  the  eclipse)  the  performance  of  the 
different  detection  algorithms  can  be  ascertained  in  a  controlled  environment.  Also, 
because  the  object  is  in  geosynchronous  orbit,  it  stays  relatively  stationary  in  the  sky,  thus 
the  object  requires  practically  no  tracking  motion  from  the  telescope  motors.  With  the 
object  location  relatively  fixed,  different  detection  algorithms  are  tested  using  the 
observations  of  the  dimming  satellite.  The  detection  algorithm  that  successfully  reports 
the  presence  most  consistently  through  the  eclipse  period  represents  the  superior 
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algorithm.  This  method  of  testing  detection  algorithms  is  far  superior  to  perfonning 
algorithm  tests  against  unknown  objects  that  may  or  may  not  actually  exist,  which 
precludes  the  possibility  of  making  a  firm  conclusion  as  to  which  algorithm  is  actually 
detecting  an  object  with  a  higher  success  rate  versus  a  detector  that  produces  more  false 
detections. 

Two  types  of  BHTs  are  compared  with  the  MHT,  one  is  the  baseline  point 
detector  used  by  the  SST  and  the  other  is  a  matched  filter  technique  (i.e.  correlation 
detector)  similar  to  that  used  by  the  Pan-STARRS  program  [41].  The  advantage  of  a 
MHT  is  gained  in  part  by  mitigating  the  aliasing  caused  by  the  undersampled  SST  images 
[43].  The  sufficient  statistics  for  both  the  BHT  and  MHT  are  derived  in  terms  of  signal- 
to-noise  ratio  (SNR)  [25,  38],  The  hypothesis  that  maximizes  SNR  while  simultaneously 
increasing  the  probability  of  detection  (Pd)  is  chosen,  thereby  providing  sub-pixel 
position  infonnation  on  the  image  location  and  increasing  Pd  over  the  BHTs. 

The  comparison  of  the  different  hypothesis  testing  methods  on  the  basis  of 
probability  of  detection  and  processing  requirements  is  made  using  data  collected  from 
the  experiment  described  in  the  next  section.  A  modeled  PSF  is  generated  using  a  phase 
retrieval  technique  presented  in  Section  4.4  and  then  it  is  dithered  to  create  the  MHT. 

This  PSF  model  is  utilized  because  it  has  been  proven  and  documented  to  work  with  the 
SST  in  the  past  [17].  In  Section  4.6  a  comparison  of  a  MHT  to  the  BHTs  is  conducted  to 
illustrate  the  advantages  of  the  MHT  as  well  as  its  additional  computational  burden.  Then 
at  the  end  of  Section  4.6,  a  derivation  shows  that  the  SNR  results  of  the  MHT  are  linearly 
related  to  the  LS  estimates  of  point  source  irradiance  which  improves  photometry. 
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4.2  The  SST  Experimental  Description 

As  stated  in  the  introduction,  the  purpose  of  this  experiment  is  not  to  find  a  new 
object  in  space  using  the  SST,  but  to  instead  run  different  detection  algorithms  on  data 
containing  a  very  dim,  but  known  object  so  that  the  relative  performance  of  different 
detection  algorithms  can  be  compared  in  a  controlled  environment.  The  experimental 
method  provides  a  data  set  that  is  used  to  form  a  clearly  supportable  conclusion  as  to 
what  algorithm  should  be  used  to  help  detect  dim  objects  in  space.  In  essence,  you 
cannot  measure  the  probability  of  detection  for  a  system  without  knowing  with  certainty 
that  an  object  is  present  to  detect.  In  addition  to  the  experimental  description,  a  basic 
overview  of  the  SST’s  design  and  current  detection  strategy  are  covered. 

4.2.1  Experimental  Setup  and  Process  Overview 

The  experiment  was  conducted  by  imaging  a  GEO  communications  satellite, 
ANIK-F1,  with  the  SST  in  a  test  mode  as  the  satellite  went  into  and  out  of  eclipse  during 
the  2012  vernal  equinox  as  illustrated  in  Figure  19.  It  was  important  to  conduct  the 
experiment  near  the  equinox  because  GEO  satellites  only  eclipse  during  that  period  of  the 
solar  cycle. 
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Figure  19.  Eclipse  experiment  overview  [41,  42]. 

There  are  many  cataloged  dim  astronomical  objects  which  could  be  used  to 
compare  detection  algorithm  performance,  but  by  imaging  eclipsing  GEO  satellites,  the 
experimental  observations  capture  both  the  effect  of  the  irradiance  division  across  pixels 
that  arises  from  objects  moving  across  the  FOV  of  the  telescope  and  the  decreasing 
irradiance  levels  of  the  satellites  as  they  enter  into  eclipse.  The  irradiance  levels  decrease 
as  the  satellites  move  through  the  penumbra  and  into  the  umbra  as  illustrated  by  the  light 
curve  produced  from  data  on  ANIK-F1  using  the  U.S  Naval  Observatory’s  (USNO)  1  m 
telescope  and  plotted  in  the  lower  right  hand  corner  of  Figure  19.  The  roll  off  of  ANIK- 
Fl’s  irradiance  during  eclipse  was  first  documented  in  a  series  of  experiments  conducted 
by  USNO  to  record  the  glint  of  GEO  satellites  shortly  before  eclipse  [50]. 

The  first  stage  in  the  experiment  is  the  collection  of  the  raw  data.  On  each  night 
of  the  experiment,  images  of  the  night  sky  containing  ANIK-F1  were  collected  using 
100ms  exposures  at  a  rate  of  8  frames  per  second.  The  telescope  was  pointed  so  that 
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ANIK-F1  was  centered  in  the  field-of-view  near  the  time  right  before  the  eclipse.  The 
orbital  elements  of  the  satellite  were  entered  into  the  SST’s  tracking  system,  so  that  the 
SST  could  be  programmed  to  follow  the  satellite  through  the  eclipse.  This  required  very 
little  tracking  movement  from  the  motors  as  the  object  being  in  GEO  orbit  appears  to  be 
stationary  in  the  sky  at  approximately  the  same  position  throughout  the  data  collection. 

As  predicted  on  many  nights  of  the  eclipse,  the  satellite  became  too  dim  to  detect  with  the 
SST’s  existing  detection  software  (described  in  the  next  section).  Once  the  data  was 
collected  it  was  recorded  and  provided  to  the  algorithm  test  team  at  AFIT  for  post¬ 
processing. 

The  next  step  in  the  experiment  was  the  pre-processing  phase.  The  raw  6144  by 
4096  pixel  SST  image  data  was  reduced  to  a  more  manageable  data  set  involving  only 
200  by  200  pixels  around  ANIK-F1.  This  allowed  for  more  efficient  use  of  memory 
resources  within  the  computer,  while  providing  a  sufficiently  large  field  of  view  to  be 
certain  the  satellite  was  fully  contained  in  the  reduced  frame  as  well  as  capturing  nearby 
stars  for  use  in  determining  the  system  PSF.  At  the  beginning  of  the  test  ANIK-F1  is 
bright  (roughly  a  magnitude  9  object)  and  clearly  visible  in  the  center  of  the  field  of  view 
of  the  telescope.  Efforts  to  manually  identify  its  position  are  further  aided  by  the  fact  that 
it  doesn’t  change  position  appreciably  throughout  the  test.  Also,  the  satellite  is  readily 
identifiable  when  it  emerges  from  the  eclipse  (again  returning  to  its  pre-eclipse 
magnitude),  thus  a  linear  trajectory  of  the  object  can  be  predicted  through  the  eclipse  and 
its  exact  position  (to  within  a  pixel)  can  be  predicted  for  every  frame,  thus  no  other  image 
registration  algorithm  is  required.  A  priori  sub-pixel  location  information  is  not  required 
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to  perform  the  experiment  as  all  existing  tests  are  designed  to  make  a  simple  binary 
decision  of  whether  the  object  is  within  the  pixel  or  not. 

The  next  step  in  the  experiment  is  to  extract  the  point  spread  function  from  the 
images  for  use  by  the  different  detectors.  Three  different  detectors  are  used  to  process 
the  data  from  this  test  for  comparison.  The  first  is  the  point  detector  (described  in  the 
next  section),  which  is  currently  used  by  the  SST,  LINEAR,  Pan-STARRS  and  other 
deep  space  object  detection  programs  [26,  41].  This  detector  does  not  utilize  a  point 
spread  function  since  it  analyzes  the  data  just  within  a  single  pixel  to  make  detection 
decisions.  The  second  detector  is  the  correlator  or  matched  filter  detector.  This  detector 
is  used  optionally  by  the  Pan-STARRS  program  to  make  detection  decisions  and  requires 
the  use  of  a  PSF  [41].  As  shown  in  Figure  20,  a  star  is  selected  to  provide  the  PSF  shape 
for  the  correlator  on  each  night.  The  selected  star  is  chosen  to  match  the  shape  of  the 
satellite  observed  near  the  start  of  the  test  in  order  to  help  maximize  the  perfonnance  of 
the  correlator.  The  correlator  as  implemented  by  the  Pan-STARRS  program  is  not 
designed  to  consider  undersampling  or  sub-pixel  motion,  so  a  single  empirically 
measured  PSF  is  used  each  night  to  implement  this  particular  detector. 


(a)  (b)  (c) 


Figure  20.  Images  of  stars  used  for  correlator  on  (a)  2012  March  13,  (b)  2012  March  14,  (c)  2012 
March  15. 
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The  PSF  model  used  for  the  MHT  requires  a  properly  sampled  PSF  that  can  be 
used  to  generate  the  PSF  shape  for  the  nine  different  hypotheses  used  in  the  test 
corresponding  to  the  nine  different  sub-pixel  locations.  This  was  done  using  a  model- 
based  approach.  This  same  optical  model  was  used  by  the  MIT  Lincoln  labs  team  to 
measure  the  PSF  in  order  to  achieve  focus  and  alignment  of  the  telescope  in  March  of 
2012  just  preceding  the  eclipse  event.  The  same  optical  model  is  used  for  characterizing 
the  PSF  in  hopes  of  leveraging  their  experience  with  the  telescope  to  compute  the  PSF 
[4].  Although  other  methods  for  computing  a  properly  sampled  PSF  from  undersampled 
imagery  exist,  performing  a  comparison  study  between  these  methods  was  not  the 
purpose  of  this  paper  [51,  52,  53,  54,  55].  Clearly,  a  better  PSF  estimate  would  lead  to 
even  better  performance  for  the  MHT  method  since  it  is  the  only  detector  tested  in  this 
study  that  utilized  of  a  properly  sampled  PSF.  The  detailed  steps  on  how  the  modeled 
PSF  used  to  construct  the  MHT  is  computed  are  described  in  Section  4.6. 

The  final  step  in  the  test  is  to  provide  each  detector  with  the  raw  image  data  in  a 
19  by  19  window  centered  on  the  pixel  containing  the  satellite  for  all  frames  of  data 
gathered  by  the  SST  of  the  satellite.  Each  detector  reports  an  SNR  for  the  satellite  for 
each  frame  of  data.  The  SNR  values  over  10  frames  are  locally  averaged  to  reduce  the 
effect  of  noise.  The  averaged  SNR  is  then  converted  to  a  probability  of  detection  for  the 
point  detector,  the  correlator  and  the  MHT  via  the  Eqs.  (2.34),  (2.35),  and  (4.20) 
respectively.  Although  a  detection  decision  could  be  made  based  on  the  reported  SNR 
for  each  detector  in  each  frame,  the  computed  probability  of  detection  reflects  the 
statistical  chance  of  making  a  correct  detection  decision  for  the  object  based  on  the 
average  SNR  and  the  estimated  noise  level  present  in  the  data.  The  computed  probability 
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of  detection  is  a  superior  perfonnance  metric  to  SNR  or  empirical  detection  frequency 
because  it  conveys  the  improvement  of  one  detector  over  another  in  tenns  that  can  be 
more  readily  understood.  Results  of  the  reported  probability  of  detection  for  each 
detector  on  each  night  are  reported  in  Section  4.7  of  this  chapter. 

4.2.2  The  SST System 

The  SST  has  a  Mersenne-Schmidt  design  selected  for  both  its  wide  FOV  and 
compactness  [13].  The  3.5  m  diameter  primary  mirror  was  built  to  meet  the  requirement 
of  detecting  small  faint  objects  with  relatively  short  integration  times,  thereby  avoiding 
streaking  of  the  satellite  image  across  multiple  CCD  pixels  so  that  the  objects  are  suitably 
modeled  as  point  sources.  Another  characteristic  of  the  Mersenne-Schmidt  design  is  a 
curved  focal  surface,  which  allows  the  SST  to  better  optimize  spot  size  across  the  field  of 
view  and  spectral  response  of  the  CCD.  Consequently,  the  curved  CCD  imager  and 
mosaic  camera  were  developed  specifically  for  the  telescope  [3]. 

4.2.3  The  SST  Detection  Process 

As  discussed  in  the  Section  2.4,  the  SST’s  baseline  detection  method  is  based  on 
the  point  detection  BFIT  used  for  Lincoln  Near  Earth  Asteroid  Research  (LINEAR) 
conducted  at  the  Experimental  Test  Site  near  Socorro,  NM.  In  nonnal  operating  mode 
the  CCD’s  15  pm  pixels  are  2  by  2  binned  and  the  array  has  a  6144  by  4096  binned 
format.  While  binning  the  data  increases  the  SNR  and  improves  detection  performance, 
it  also  increases  the  amount  of  undersampling  of  the  data  [43]. 
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4.3  Binary  Hypothesis  Testing  (BHT) 

Replacing  the  baseline  BHT  with  a  correlation  based  BHT  improves  the  SST’s 
detection  perfonnance.  This  method  for  improving  the  SST’s  detection  performance  is 
explored  by  [16]  in  which  a  comparison  of  the  SNR  of  two  different  binary  hypothesis 
tests,  a  correlator  and  a  point  detector,  is  made  using  the  ANIK-F1  experimental  data.  In 
this  comparison,  stars  in  the  FOV  are  used  to  estimate  the  total  system  PSF.  Figure  20 
(a-c)  shows  irradiance  maps  of  each  star  within  a  19  by  19  window  cropped  from  the  SST 
images  on  three  consecutive  nights  and  used  in  the  correlation  detector.  One  item  to  note 
is  that  in  Figure  20  (b  &  c)  the  star  appears  to  have  a  different  shape  than  the  star  in 
Figure  20  (a).  The  change  in  apparent  shape  of  the  star  is  due  to  the  images  being 
centered  at  different  sub-pixel  locations.  If  the  shape  of  the  object  of  interest  is  not  the 
same  as  the  total  system  PSF  used  in  the  correlator,  the  detection  perfonnance  will  be 
degraded. 

The  SST’s  threshold  used  for  detection  during  the  technical  demonstration  period 
is  y  =  6,  which  inherently  sets  the  probability  of  false  alann.  The  probability  of  false 
alann  is  defined  as  the  chance  that  a  pixel  that  contains  only  background  light  (no  object) 
will  produce  a  detector  output  that  exceeds  the  threshold  value  of  6.  When  objects  are 
not  present  in  the  pixel  the  operations  described  in  Eq.  (2.34)  and  Eq.  (2.35)  are  designed 
to  produce  unit  variance  zero  mean  Gaussian  random  variables.  Therefore,  the 
probability  of  false  alanns,  PFA,  per  pixel  is 

PFA=P{SNRBaseline>6\H0) 

=  P(SNRcon>6\H0)  (4.1) 

e  2  dt  =  9.87e  -  010, 
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where  SNRBaseline  is  the  output  of  the  point  detector,  SNRcorr  is  the  output  of  the  correlator, 
and  H0  is  the  hypothesis  that  no  object  is  present  in  the  pixel. 

4.4  The  SST  PSF  Modeling 

The  long  exposure  modeled  PSF,  hL  (x)  centered  on  a  pixel  is  then  computed 
from  Eq.  (2. 14)  where  x  is  the  Nyquist  pixel  coordinates.  One  important  property  of  the 
Nyquist  sampled  PSF  is  that  a  sub-pixel  shift,  Ax,  of  the  model  does  not  significantly 
change  its  shape.  Modeling  the  effects  of  Ax  on  the  image  irradiance  pattern  is 
necessary  because  the  point  source  is  not  always  in  the  center  of  a  pixel.  To  reproduce 
the  change  in  irradiance  pattern  measured  by  each  of  the  30  pm  pixels  in  the  CCD  as  a 

function  of  Ax,  the  modeled  PSF  is  down-sampled  using  the  ratio,  g,  between  the  30  pm 
pixels  and  the  Nyquist  pixels  size  from  Eq.  (2.38).  The  shifted  and  down-sampled  PSF  is 

hsamp  (m,  Ax)  -  j  h,  (x)  •  8 (gm  -  x  -  gAx) dx,  (4.2) 

thus  the  sampled  irradiance  is 

Kamp Of  tsx)  =  0-  hsamp (m.  Ax)  +  B,  (4.3) 

where  B  is  the  background  light,  and  6  is  the  total  number  photons  emitted  from  the 
object  per  integration,  and  m  is  a  integer  valued  pixel  location  in  the  CCD  array.  To 
build  the  PSF  model,  estimates  of  the  coefficients  Z5  -  Zn  were  made  using  the  Donut 
algorithm  and  inserted  into  the  PSF  model  using  Eqs.  #  (2.3)-(2.6)  [17]. 

The  method  of  LS,  described  in  Chapter  III  using  Eq.  (3.37),  is  used  to  jointly 
find  Z4  and  r0  from  a  star  selected  the  first  frame  in  the  ANIK-F 1  experimental  data 
(see  Figure  20).  A  color  map  of  the  PSF  model  generated  using  a  the  SST  star  image  on 
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2012  March  14  with  both  Nyquist  and  30  jam  sampling  are  shown  respectively  in  Figure 
21  (a)  and  (b).  In  the  SST  data,  d(m),  the  star  image  can  be  centered  at  any  sub-pixel 
location  and  the  corresponding  irradiance  pattern  changes.  By  shifting  the  modeled  PSF 
the  changes  in  the  irradiance  pattern  in  the  star  data  can  be  captured.  If  the  PSF  is 
undersampled  spatially,  as  is  the  case  with  the  30  pm  detected  PSF,  the  shifted  PSF  will 
have  aliasing  artifacts  [43].  Figure  21  (c)  depicts  the  aliasing  artifacts  produced  when  the 
undersampled  PSF  is  shifted  using  the  Fourier  transform  shift  method 

v(m)  =  cos(m-Ax)  +  zsin(m-  Ax),  and  (4.4) 

Khift  (m,Ax)  -  RE  { J  {hsamp  (m)[  •  u  (w)}] ,  (4.5) 

where  Ax  magnitude  of  sub-pixel  PSF  shift  such  that  the  modeled  irradiance  is 

Khift  {m,Ax)  =  6  •  hshift  (m.  Ax)  +  B.  (4.6) 

However,  if  the  Nyquist  sampled  model  is  shifted  before  down  sampling  to  30  pm  using 
Eq.  (4.2)  as  shown  in  Figure  21  (d)  and  (e)  the  irradiance  pattern  does  not  have  the  same 
aliasing  artifacts  as  seen  in  Figure  21  (c). 
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Figure  21.  The  SST’s  phase  retrieved  PSF  on  2012  March  14.  (a)  Centered  Nyquist  sampled  PSF, 
a  =  2.75e-07  m,  on  a  2048  by  2048  grid  (b)  Centered  down-sampled  PSF,  a  =  30  pm,  on  a  19  by 
19  grid  (c)  Model  PSF,  a  =  30  pm,  shifted  to  the  lower  right  hand  comer  of  pixel  (10,10)  with 
aliasing  artifacts,  (d)  Nyquist  sampled  model  PSF  shifted  without  aliasing  artifacts,  (e)  Sampled 
model  PSF,  a  =  30  pm,  up  sampled  from  Nyquist  sampled  model  PSF  shifted  to  the  lower  right 
hand  comer  of  pixel  (10,10)  without  aliasing  artifacts. 


The  accuracy  of  the  sampled  irradiance  models  shifted  two  different  ways  can  be 
quantified  using  the  correlation  coefficient  [56].  The  correlation  coefficient  measures 
how  accurately  the  modeled  irradiance  pattern  matches  the  measured  irradiance  pattern 
on  a  scale  from  zero  to  one,  where  a  value  of  one  means  they  are  perfectly  correlated  and 
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zero  means  they  are  uncorrelated.  The  correlation  coefficient  between  both  irradiance 
models,  Eq.  (4.6)  and  Eq.  (4.3),  and  the  measured  data  are  computed  respectively  as 


and 


P shi/i  (M  =  p{d  (m),ishift(m,Ax)} 


Mj 

I 

m= 1 

(<f  (m)-B'j 

(e-hshift(mM)) 

^ZidN-B)2^ 

\Md 

Ax 

m- 1 

))! 

Psamp  (Ax)  =  p{d  (m),isamp(m,  Ax)} 

Mj 

Y\(d{m)-B)(d-  hsamp  (m,  Ax)) 

_  _ m= 1 _ 

Md  (M~d  ~ 

^0-hsamp(m, Ax)) 

\  m= 1  \  m= 1 


(4.8) 


Figure  22  is  a  plot  of  the  maximum  values  of  pshift  (Ax)  and  p  (Ax)  for  images 

of  ANIK-F1  on  2012  March  14  as  its  irradiance  is  split  between  pixels.  As  ANIK-F1  ’s 
irradiance  moves  between  pixels,  even  though  ishift(m,Ax)  is  shifted  to  maximize  the 

correlation  coefficient,  the  pshift{ Ax)  goes  down.  In  contrast,  the  maximum  value  of 
p  (Ax)  is  relatively  constant  regardless  of  where  in  the  pixel  the  irradiance  of  ANIK- 

F 1  is  located,  thus  illustrating  the  importance  generating  a  Nyquist  sampled  model.  In 
addition,  the  strong  correlation  between  the  modeled  PSF  and  the  data  indicates  that  the 
model  is  an  accurate  representation  of  the  SST  PSF. 
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Figure  22.  Correlation  between  two  different  irradiance  models  and  the  images  of  AN1K-F1  on 
2012  March  14  as  its  irradiance  is  split  between  pixels.  One  model  is  shifted  on  the 
undersampled  grid,  pMfi  (Ax),  and  the  other  on  a  Nyquist  grid  and  then  down  sampled, 

Psamp  (4*). 


4.5  Data  Normalization  Using  Outlier  Rejection  Techniques 

Another  feature  of  the  proposed  algorithm  design  is  that  background  noise 
statistics  are  computed  using  a  reduced  set  of  data  from  the  window  around  the  pixel  to 
be  tested.  This  new  noise  power  estimation  technique  has  the  feature  that  it  rejects  any 
noise  sample  in  the  window  surrounding  the  pixel  to  be  tested,  whose  values  do  not 
conform  to  those  predicted  by  Gaussian  statistics.  In  this  way,  bad  pixels  and  nearby 
stars  are  not  used  to  compute  the  noise  generated  by  the  background  light.  Current 
algorithms  used  by  the  SST  and  LINEAR  use  all  the  pixels  in  the  window  surrounding 
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the  point  to  be  tested  to  compute  the  local  noise  standard  deviation,  a  [25].  In  this 
process  the  background  B  is  computed  as  in  Eq.  (2.23).  The  squared  deviations,  /,  from 
the  background  within  the  window  are  computed  as 

X(m)  =  (d(m)-B)2 .  (4.9) 

These  squared  deviations  follow  a  chi-squared  distribution  based  on  the  assumed 
Gaussian  nature  of  the  data,  d.  Since  samples  size  is  large  (Mj  =  36 1  j  the  distribution  is 

symmetric.  Therefore,  the  mean  of  the  squared  deviations,  M,  within  the  window  is 
found 


M  =  E\_x(m)~\, 


(4.10) 


then  the  standard  deviation,  S,  of/  is  computed: 


5 


Md 


YjX2(™) 


m=\ 


-- M 


(4.11) 


A  new  noise  standard  deviation,  is  computed  from  the  window  using  Eq. 
(2.24)  by  excluding  any  pixel,  m,  in  the  calculation  where  %  ( m )  >  (M  +  3  ■  S).  The  new 
noise  standard  deviation  is  included  in  the  following  MHT  for  improved  detection 
performance  by  normalizing  the  data  as  in  Eq.  (2.34)  by  replacing  crwith  d- 


4.6  Multi-hypothesis  testing  (MHT) 

A  multi-hypothesis  (M-ary)  detector  is  introduced  because  the  image  of  a  space 
object  does  not  always  fall  in  the  center  of  the  pixel.  In  addition,  simple  correlation 
operations  are  not  desirable  because  the  shape  of  the  sampled  PSF  changes  depending  on 
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where  the  object  is  imaged  on  the  array.  In  order  to  account  for  the  possibility  that  the 
image  is  in  different  places  within  the  detector,  we  introduce  a  multi-hypothesis  test 
(MHT)  strategy.  The  hypothesis  that  an  image  of  a  space  object  is  not  present  in  the 
pixel,  H0,  plus  the  nine  different  sampled  PSF  shapes,  shown  in  Figure  23  fonn  the  ten 
hypotheses  for  the  MHT,  This  choice  of  hypotheses  captures  a  great 

deal  of  the  spatial  dependence  of  the  PSF  while  only  introducing  one  order  of  magnitude 
more  computations. 


Figure  23.  Hypothesis  that  the  point  source  image  is  in  either  the  center  of  the  pixel,  Hi, 
on  the  sides,  H2-H5i  or  corner  of  a  pixel,  H6-H9. 
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According  to  Kay,  the  M-ary  decision  to  select  one  hypothesis,  Hk,  over  another 
hypothesis,  Hp  is  the  M-ary  maximum  likelihood  (ML)  decision  rule.  Kay  derives  the 
ML  decision  rule  using  a  Bayes  risk,  R ,  approach  under  the  assumption  of  uniform  cost 
and  equal  priors  to  decided  between  the  different  hypotheses  He  assigns 

cost,  Cik  to  the  decision  if  Hj  is  selected  when  Hk  is  true.  The  risk  is  calculated  as 


R  =  YLC>Ah< \Ht)p(H,),  (4.12) 

i= 0  j= 0 

where 


JO  /  =  k 
{ 1  i  k 


(4.13) 


By  assuming  uniform  cost,  each  of  the  hypotheses  is  given  the  same  priority.  With  this 
choice  of  cost,  the  information  provided  on  position  of  the  object  in  the  pixel  is  weighted 
the  same  as  the  binary  decision  of  whether  an  object  is  detected  or  not  in  the  MHT.  The 
position  information  from  the  MHT  is  important  for  improving  the  SST’s  metric  accuracy 
over  the  BHT.  However,  detection  is  more  important  than  position  accuracy  and  this 
choice  of  cost  could  reduce  the  MHT’s  detection  performance.  Even  with  this  penalty  in 
detection  perfonnance,  the  MHT  will  be  shown  to  still  outperform  the  correlator  and  the 
baseline  detector  in  tenns  of  probability  of  detection.  In  addition,  by  making  this  choice 
of  cost,  the  derivation  of  the  detector's  sufficient  statistic  is  greatly  simplified  and  likely 
produces  a  detector  that  is  more  computationally  efficient  to  apply. 

Kay  shows  that  to  minimize  R,  the  hypothesis  that  minimizes  the  average  cost  of 

deciding  Hj  if  d  (w,z)  is  observed,  C.  {d  (w,z)),  should  be  selected  where 
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(4.14) 


Ci(d(w,z))  =  YJCikp(Hk  \d(w,z)), 

k= 0 

for  i  =  0,1, ....,9 .  Inserting  Eq.  (4.13)  into  Eq.  (4.14) 


Ci(d(w,z))  =  YJp(Hk\d(w,z)) 

k= 0 
k^i 

=  lLP(Hk  I d(w’z))~p{Hi  I d(w,z)). 

k= 0 


(4.15) 


Since  p  ( H k  \  d  (w,  z))  is  not  a  function  of  i  the  risk  is  minimized  by  maximizing 
p{Hi  |  d  (w,z)).  Therefore  the  minimum  risk  decision  rule  is  choose  Hk  if, 


p{Hk  |  c/(w,z))  >  p{Ht  |  c/(w,z))  V/  ^  k. 


(4.16) 


Then  for  equal  priors  probabilities, 


p[Hi\d(w,z)) 


_  p(d(w,z)\Hi)p(Hi ) 
p(d(w,z )) 


/7(flf(w,z)|i7f) 


10 


p(d(w,z )) 


(4.17) 


so  maximizing  p[d{w,z)  \  //,.)  maximizes  p(//(  |  <7(rv,z))  and  the  ML  decision  rule 


becomes  choose  Hk  if  [38] 


4  p(^(w,z)|if/r)  p(y(w,z)|//,.) 

Gt  p{d(w,z)  \H0)  p(d(w,z)\H0) 


Assuming  the  prior  probability  that  an  object  is  not  in  a  pixel  is  equal  likely  as  an  object 
being  in  the  pixel  is  not  a  precise  choice  due  to  the  density  of  stellar  objects.  However, 
the  choice  of  equal  priors  is  chosen  because  the  true  probabilities  are  unknown. 
Furthennore,  the  baseline  detector  and  correlator  make  this  same  assumption  so  the 
comparison  of  the  MHT  to  these  detectors  is  a  like  comparison.  As  with  cost  selection,  a 
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better  choice  in  priors  could  result  in  a  better  perfonning  detector,  but  choices  of  uniform 
cost  and  equal  priors  are  shown  to  produce  a  MHT  that  outperforms  the  existing  BHT. 

Since  Eq.  (4. 18)  has  the  same  forms  as  the  LRT  given  in  Eq.  (2.25),  the  sufficient 
statistic  for  BHT  given  in  Eq.  (2.34)  can  be  applied  to  the  MHT  accounting  for  the 
additional  hypotheses  based  on  sub-pixel  position  shifts  listed  in  Table  1,  at  &  /?.. 
Therefore,  the  ML  decision  rule  can  be  determined  as  a  function  of  SNR  for  each 
location,  /  (57W?.),  so  that  Eq.  (4.18)  becomes  chose  Hk  if 

A„i=f(SNRl)>f(SNRl)=A(li  Vi*k.  (4.19) 

where  MHT  sufficient  statistic  is 


M,  M, 


E  Z  (rf(  W’ Z)  -  R)hsamp  ( W  +  Cx~ai,Z  +  Cy  ~  fi. )  H‘ 


SNR .  =  ■ 


M,,  M,, 


]hsamp2(W’Z) 


> 

< 


I'm- 


ary 


(4.20) 


•  H  :  Hypothesis  that  no  satellite  is  present. 

•  H :  Hypothesis  that  a  satellite  is  present  (see  Table  2) 


Table  2.  Alternative  Hypothesis  Sub-pixel  Shifts  (corresponding  to  Figure  23) 


Alternative  (i) 

Horizontal  Shift  (a;) 

Vertical  Shift(Pi) 

1 

0 

0 

2 

0 

-15pm 

3 

0 

15  pm 

4 

15  pm 

0 

5 

-15pm 

0 

6 

15  pm 

-15pm 

7 

15  pm 

15  pm 

8 

-15pm 

-15pm 

9 

-15pm 

15  pm 

16 


Applying  a  generalized  likelihood  ratio  test  (GLRT),  PSF  shift  estimates, j , 
used  to  determine  if  there  is  a  detection  in  the  ML  location  [38] 

p(d(w,z)\{Hl,H2,...,H9})> 


are 


p(d(w,z)\H0) 

p[d(w,z)  | 

p(d(w,z)\H0) 


< 

H0 

Ht 

> 

< 

Hn 


r. 


where 


( # )  =  arg max  ( d  ( w,z)\  7/ ) . 

V  ’  a,  ,#,1=1:9 


Therefore,  the  GLRT  can  be  rewritten  as 

:p{d(w,z)\Hi) 
p(d(w,z)\H0 ) 


max , 

i=l:9 


H, 

> 

< 


7 


=  max 

i=l:9 


p(d(w,z)\Hi) 

p(d(w,z)\H0) 


7- 


(4.21) 


(4.22) 


(4.23) 


With  this  approach,  the  M-ary  hypothesis,  HM_arv,  that  satisfies  the  ML  decision  rule  and 
and  exceeds  the  detection  threshold,  yM_m  ,  is  detennined  by  finding 

TJ 

n  M-ary 

max  {SN II  )  ^  yM_aty,  (4.24) 

H0 

which  provides  sub-pixel  image  location  infonnation  and  detection  simultaneously. 

An  important  goal  in  deriving  the  MHT  is  to  improve  the  probability  of  detection, 
PD,  without  raising  the  probability  of  false  alarm,  PFA .  In  this  case  the  probability  of  false 
alann  is  the  chance  that  a  pixel  with  no  space  object  in  it  will  be  classified  as  having  one. 
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In  mathematical  terms,  this  is  the  probability  that  the  SNR  output  of  the  detector  with  no 
object  present  will  exceed  the  detection  threshold.  Computing  the  false  alarm  probability 
is  a  challenging  task  and  may  prove  mathematically  intractable.  Instead  we  compute  an 
upper  bound  on  the  probability  of  false  alann  and  then  use  that  upper  bound  as  our 
estimate  of  the  false  alarm  probability.  This  guarantees  that  the  new  MHT  test  will  not 
raise  the  probability  of  false  alarm  over  the  existing  BHT  used  by  the  SST  and  LINEAR. 

Extending  the  MHT,  from  testing  only  one  pixel  for  H0  -  H9  to  testing  each  pixel 

in  the  frame,  results  in  repeating  the  same  test  four  times  in  the  corners  and  two  times  on 
the  sides  of  each  pixel  as  depicted  in  Figure  24.  Therefore,  to  minimize  the  number  of 
hypothesis  tests  each  pixel  only  needs  to  be  tested  in  the  center,  on  one  corner  and  on  two 
sides  and  the  overlap  of  the  grids  forms  the  9  hypotheses  for  each  pixels.  While  this 
technique  reduces  the  MHT’s  number  of  computations,  it  is  a  slight  departure  from 
testing  each  pixel  9  times  because  the  window  statistics  around  each  pixel  varies. 


Figure  24.  Illustration  depicts  the  overlap  of  the  comers  and  sides  of  the  pixels  on  the  CCD.  (The 
black  dots  represent  the  center  of  the  pixel;  the  green  dot  is  the  comer;  blue  and  red  dots  are  the 
sides) 
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Minimizing  the  number  of  hypothesis  tests,  thereby  decreasing  processing  time  and  PF4. 

The  MHT  computational  cost  computed  in  terms  of  the  number  of  floating  point 
operations  (Flops)  is  ~6  times  more  than  the  baseline  approach  as  listed  in  Table  3. 


Table  3.  Summary  of  the  Flops  required  for  Each  Detector 


Parameter 

Baseline  Flops 

Correlator  Flops 

M-ary  Flops 

B 

Md2 

Md2 

4  Mi 

a 

3  Mi 

3  Mi 

12  Mi 

SNR 

2 

2Mi+2 

mA-+2 

Total  Flops 

4Ma2+2 

6Mi+2 

24Mi+2 

Two  simplifying  assumptions  are  made  to  find  the  upper  bound  of  the  PFA  for  the 
M-ary  test.  The  first  is  considering  PFA  for  each  alternative  hypothesis  of  the  M-ary  test 

4 

to  be  mutual  exclusive  such  that  Pi  Ha  =  0.  The  second  assumption  is  that  the  result  of 

<2=1 


each  individual  alternative  in  M-ary  test  is  statistically  independent  of  each  other.  Under 
those  two  conditions  the  PFA  can  be  bounded  above  by  extending  Eq.  (4. 1)  to 


p  =  p 

1  FA  1 


<P 


f  4 

u 

V  <2=1 
f  4 


f  i 

\H0 

-P 

fRI^o 

J 

J 

(4.25) 


=  I H0)  =  4 .  P(SVKm_„  >  6 1  H„) 


a= 1 


( 4  x  9.87e  -  010  =  3.94e-009. 


79 


The  estimated  PFA  is  higher  for  the  M-ary  test  than  for  the  BHT,  but  it  can  be  reduced  by 
raising  the  M-ary  detection  threshold  to  yM_ary  =  6.2212  so  that  PFA  «  4  x  2.467e-010 
*  9.87e-010. 

An  example  of  the  detection  perfonnance  gains  from  the  M-ary  test  are  shown  in 
Figure  25(a-f)  using  data  from  the  satellite  eclipse  experiment.  To  produce  the  plots, 
running  averages  with  a  50  frame  window  for  the  baseline  detector  SNR,  juBaseline,  the 

correlator  SNR,  jucorr,  and  M-ary  test  SNR,  nM  ,  were  found  in  the  threshold  region 

for  all  three  detectors.  The  probability  of  detection  for  the  baseline  detector,  the 
correlator,  and  M-ary  test  as  a  function  of  running  average  are  respectively 

^  1  -(t-MBaseline)2  'j 

PDBaseKn\P Baseline)  =  J  HTZ  ^  ^  ’  (4-26) 

V  6  V2^-  ) 

1  ~( t-Mcorrf  V 

(Peon  )  =  J  ~nr=e  2  dt  ,  and  (4.27) 

v6  v2^"  y 

/  2  \  3 

PDM_ary{PM-ary)~  j  ^  ’  (4.28) 

Kb.22\2^Ln  j 

based  on  the  Gaussian  noise  assumption  and  the  fact  that  tracklets  require  three 
consecutive  frames  for  detection  [25],  Note  that  the  threshold  in  Eq.  (4.28)  is  adjusted  to 
keep  the  PFA  approximately  the  same  for  all  three  detectors. 

The  probability  of  detection  computed  using  Eqs.  #  (4.21),  (4.22)  and  (4.23)  are 
based  on  the  running  average  SNR  output  from  each  detector  and  represent  the  average 
probability  of  detection  that  can  be  expected  as  the  satellites  irradiance  decreases  based 
on  the  Gaussian  distribution  of  the  noise.  In  the  case  of  the  M-ary  test,  the  running 
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average  is  produced  using  the  highest  instantaneous  SNR  output  from  each  of  the  four 
alternatives  in  the  M-ary  test  per  pixel,  which  is  selected  based  on  the  ML  decision  rule 
using  Eq.  (4. 12).  In  this  sense,  in  terms  of  detection  the  results  are  interpreted  as  a  binary 
detection  solution;  either  detection  occurred  in  the  pixel  or  it  did  not.  Having  a  single 
probability  of  detection  for  each  CCD  pixel  tested  with  the  MHT  is  important  so  that  it 
can  be  compared  with  the  BHTs.  The  other  important  information  resulting  from  the  M- 
ary  test  about  the  sub-pixel  location  necessary  for  improving  the  SST’s  metric  accuracy  is 
not  included  in  the  calculation.  In  addition,  the  fact  four  tests  are  done  at  each  pixel 
producing  SNR  values  that  each  may  exceed  the  detection  threshold  is  not  factored  into 

the  PDm  ( jUM_ary )  calculation.  Including  the  additional  test  results  in  the  calculation 

would  only  serve  to  raise  the  probability  of  detection  of  the  M-ary  test  because  each  of 
the  additional  hypotheses  could  also  result  in  the  detection  of  the  object.  Therefore,  Eq. 
(4.23)  only  provides  the  lower  bound  estimate  of  probability  of  detection  for  the  M-ary 
test. 

When  ANIK-Fl’s  irradiance  is  high,  all  three  detectors  can  detect  the  satellite,  but 
as  the  satellite  dims  as  it  enters  eclipse  the  detector  perfonn  differently.  Eventually,  the 
satellite  becomes  so  dim  that  it  is  undetectable  by  any  of  the  detectors.  The  area  of 
interest  then  becomes  the  detection  threshold  hold  region.  As  seen  in  Figure  25,  on  all 
six  nights  the  M-ary  detector  detection  performance  significantly  exceeds  both  the 
correlator  and  the  baseline  detector.  This  means  that  the  M-ary  test  detects  much  dimmer 
objects  that  the  baseline  algorithms.  The  perfonnance  gains  seen  on  2012  Mar  23  are 
only  due  to  better  calculations  of  the  background  noise  statistics.  On  that  date  there  was 


81 


a  narrow  total  PSF  and  a  lack  of  aliasing  because  the  SST  tracked  the  object  in  the  same 
pixel. 


(a)  (b) 


(c)  (d) 


(e) 


(f) 


Figure  25.  Comparison  of  the  baseline  detector,  the  correlator,  and  M-ary  test  probability  of 
detecting  (Pd)  AN1K-F1  as  it  enters  eclipse  on  (a)  2012  March  13  (b)  2012  March  14  (c)  2012 
March  15  (d)  2012  March  21  (e)  2012  March  22  (f)  2012  March  23  where  the  PFA  is  equal  for  all 
three  detectors. 


Using  data  from  six  nights  of  ANIK-F1  eclipse  observations,  a  direct  comparison 
between  the  baseline  detector  and  both  the  correlator  and  M-ary  test  shows  an  improved 

probability  of  detecting  a  space  object.  A  plot  of  Pr>  (<ucorr)  and  PD  {pM-ary) 

versus  P  ( u„  ..  j  for  the  nights  of  2012  Mar  13-15  and  2012  March  21-23  is  shown 

D Baseline  \  r* Baseline  )  0 

in  Figure  26.  The  plot  shows  that  the  M-ary  test  has  a  higher  probability  of  detection  that 
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the  correlator  and  the  baseline  detector  at  the  detection  threshold.  On  2012  March  23,  the 
performance  gains  were  not  as  dramatic  because  the  objects  irradiance  was  concentrated 
in  one  pixel  due  to  excellent  tracking  of  the  satellite.  However,  when  the  object  is  not  is 
not  centered  on  a  pixel  the  M-ary  test  has  a  PD  when  the  baselines  has  only  a 

PDg  =0.5.  This  30-50%  demonstrated  improvement  in  the  probability  of  detection 

means  that  significantly  more  dim  objects  like  small  asteroids  will  be  found  with  the  SST 
using  the  M-ary  detector. 


Figure  26.  Composite  plot  of  the  probability  of  detecting  AN1K-F1  as  it  enters  with  either  the 
correlator  or  M-ary  test  versus  the  baseline  detector  for  the  nights  of  2012  Mar  13-15  and  2012 
March  21-23. 


In  addition  to  improved  detection  performance,  the  M-ary  test  can  also  provide 
better  estimates  for  object  irradiance  than  the  baseline  detector.  The  baseline  detector 
estimates  are  made  by  adding  up  the  number  of  digital  counts  in  the  pixels  where  the 
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object  was  detected  to  estimate  the  object’s  irradiance,  0Haseline.  In  contrast,  the  SNR 
output  of  the  M-ary  test  is  linearly  related  to  the  LS  estimate  of  the  objects’  irradiance, 

A 

@M-ary’  in  terms  of  digital  counts  of  by  substituting  Eq.  (4.20)  into  Eq.  (3.33) 


e 


SNR 


M-ary 


■c 


M-ary 


\Yjhsamp(m) 


(4.29) 


Figure  27  shows  that  by  using  the  M-ary  test  results  to  estimate  the  object 
irradiance  a  much  higher  digital  count  is  determined  because  it  is  counting  information 
across  the  19  by  19  window  weighted  by  the  PSF  rather  than  only  the  pixels  exceeding 
threshold.  Estimating  the  digital  count  more  accurately  should  improve  the  SST’s 
photometric  fit  results,  which  have  been  shown  to  have  a  higher  variance  for  objects  with 
magnitudes  fainter  than  18  Mv  [3]. 


Figure  27.  Comparison  of  the  baseline  detector  and  the  M-ary  test  (least  squares)  estimates  of 
AN1K-F1  irradiance  on  2012  March  14  as  it  enters  eclipse. 


84 


4.7  Conclusions 


A  MHT  can  provide  a  significant  improvement  in  the  SST’s  detection  capability 
over  a  BHT  when  the  object  is  not  centered  on  a  single  pixel,  as  shown  in  Figure  26. 
Since  the  SST’s  mission  is  to  detect  unknown  space  objects,  rather  than  track  known 
objects  it  is  unlikely  that  any  object  will  repeatedly  fall  in  the  center  of  a  pixel  during  the 
three  consecutive  frames  used  for  detection.  Recalculating  the  standard  deviation  of  the 
background  light  in  the  window  surrounding  the  pixel  being  tested  by  not  including 
outliers  also  contributes  to  the  MHT’s  improved  Pd.  The  gains  in  detection  perfonnance 
by  the  MHT  are  realized  by  mitigating  the  aliasing  effects  of  undersampling  using  phase 
retrieved  PSF  model,  but  come  at  the  cost  of  a  600%  increase  in  processing  power. 
Simultaneous  to  detection,  the  MHT  also  provides  sub-pixel  position  information  and 
more  accurate  estimates  of  object  irradiance.  These  three  improvements  to  the  SST’s 
perfonnance  by  the  MHT  come  with  a  manageable  computational  cost  that  can  be 
afforded  with  relatively  inexpensive  modem  computers  compared  to  the  cost  of  enhanced 
optics  and  therefore  give  good  cause  to  investigate  the  implementation  of  this  detector. 
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V.  Phase  Retrieval  with  a  Short  Exposure  Atmosphere 


Since  the  SST  is  considering  a  new  camera  that  can  do  short-exposure  imagery,  the 
phase  retrieval  algorithm  described  is  designed  for  the  primary  purpose  of  diagnostics  of 
an  optical  system  using  point  source  image  intensity  data  (i.e.  system  impulse  response). 
Two  cases  are  considered  first  with  no  atmosphere  present  and  the  second  with  short 
exposure  imagery  data.  The  intended  application  of  this  technique  is  to  provide 
quantitative  feedback  during  the  focus  and  alignment  of  large  three  mirror  telescopes  like 
the  Space  Surveillance  Telescope  (SST),  James  Webb  Space  Telescope  (JWST),  and 
Large  Synoptic  Survey  Telescope  (LSST)  [17,  21,  44],  For  these  telescopes  phase 
aberrations  must  be  recovered  using  post  processing  because  the  focal  surface  array 
cannot  measure  the  phase  directly. 

The  primary  advantage  of  this  phase  retrieval  technique  over  existing  phase 
retrieval  techniques  is  that  it  converges  on  the  correct  set  of  Zernike  coefficients  while 
the  telescope  is  in-focus.  For  three  mirror  telescopes  this  is  an  important  capability  since 
defocusing  the  telescope  requires  movement  of  both  the  secondary  and  tertiary  mirrors 
rather  than  only  translating  the  camera  about  the  focal  plane  or  with  the  addition  of  a  lens 
for  defocusing  [17,  46].  In  addition,  in-focus  phase  retrieval  would  enable  optical 
diagnostics  of  the  telescope  using  its  standard  imagery  data  so  the  status  of  alignment 
could  be  monitored  continually  from  star  data  in  the  field  of  view  (FOV).  The  new  phase 
retrieval  algorithm  described  in  this  chapter  achieves  in-focus  phase  retrieval  by 
estimating  the  Zernike  coefficients  using  both  an  estimated  image  plane  electric  field  (E- 
field)  and  the  measured  intensity  data.  Use  of  the  estimated  E-field  rather  than  the 
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estimated  phase  avoids  the  problems  associated  with  unwrapping  2-D  phase  produced  by 
the  estimation  algorithm. 

This  chapter  is  broken  up  into  four  sections.  The  introduction  covers  the  tools 
used  for  propagating  the  field  with  a  Fourier  transform  of  the  generalized  pupil  function, 
estimating  the  point  spread  function  (PSF)  from  point  source  data,  and  estimating  the 
electric  field  with  the  Gerchberg-Saxton  algorithm.  Then  Section  5.2  demonstrates  the 
advantages  of  estimating  Zemike  coefficients  with  the  E-field  by  comparing  the 
correlation  of  E-fields  and  the  correlation  of  intensity  patterns  produced  by  different 
phase  aberration.  Within  this  section  a  new  phase  retrieval  algorithm  and  its  subroutines 
is  described.  This  section  also  covers  simulations  are  used  to  show  the  performance 
gains  of  the  new  algorithm  with  a  wide  variety  of  phase  aberrations  compared  to  the 
intensity  based  LS  estimation  technique  used  in  isolation.  Then  Section  5.3  covers  two 
laboratory  demonstrations  that  are  used  to  show  how  the  new  phase  retrieval  algorithm 
performs  with  measured  data.  The  results  support  the  conclusion  that  the  new  algorithm 
perfonns  well  at  estimating  Zernike  coefficients  for  phase  aberration  of  an  optical  system 
in  or  near  focus. 

5,1  Introduction 

The  following  background  is  intended  to  describe  the  tools  needed  for  phase 
retrieval  including:  how  the  fields  are  propagated  with  a  Fourier  transfonn  of  the 
generalized  pupil  function,  estimating  the  point  spread  function  (PSF)  from  point  source 
data,  and  estimating  the  electric  field  with  Gerchberg-Saxton.  This  summary  of 
techniques  provides  the  fundamental  background  for  understanding  how  the  new  phase 
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retrieval  algorithm  is  implemented.  One  of  the  most  important  techniques  is  the 
propagation  of  the  field  through  a  Fourier  transform  of  the  generalized  pupil  function, 
which  is  critical  for  each  step  in  the  phase  retrieval  process.  Since  the  data  contains 
noise,  the  expectation  maximization  algorithm  described  for  PSF  estimation  from  the  data 
is  a  critical  input  to  the  new  Zemike  coefficient  estimation  technique.  The  other 
technique  critical  for  the  new  in-focus  Zemike  coefficient  phase  retrieval  algorithm  is  the 
estimation  of  the  electric  field  from  the  estimated  PSF,  which  is  done  using  Gerchberg- 
Saxton  (GS)  phase  retrieval. 

5.1.1  Generalized  Pupil  Function 

This  section  is  a  recap  of  the  description  of  the  generalized  pupil  function  from 
Chapter  II  intended  to  highlight  the  dependence  of  the  image  plane  electric  field  and 
intensity  on  the  wavefront  error,  W,  in  the  exit  pupil  of  an  optical  system.  The  majority 
of  the  wavefront  error  in  an  optical  system  can  be  decomposed  into  /V- number  of  Zemike 
polynomials,  ^  -  (f)N ,  represented  as 

W(ux)  =  Zx-</\  (ux)  +  ...  +  ZN  ■  tf>N  (w,),  (5.1) 

where  w,  is  a  coordinate  in  the  pupil  plane  and  Z,  —  ZN  are  the  Zernike  coefficients  [31]. 
The  aberrations  can  be  represented  in  the  generalized  pupil  function,  T(ux),  as 

T(ux)  =  A(u  j)e;'rW,  (5.2) 

where  A(ux)  is  the  pupil  transmittance  function.  From  the  generalized  pupil  function  a 
model  of  the  electric  field  as  a  function  of  image  plane  pixel  coordinates,  m,  and 
wavefront  error,  is  computed  using  a  discrete  Fourier  transform, 
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H(m,W)  =  YJP{ul)eJ2™“'  =  ?\_P{ux)_ 


(5.3) 


“i 

and  the  corresponding  PSF,  h(m,W),  is 

h(m,W)  =  \H(m,W)f  =H(m,W)H(m,W)* .  (5.4) 

5.1.2  Intensity  Model  and  PSF  Estimation 

To  conduct  Gerchberg-Saxton  (GS)  phase  retrieval,  an  accurate  PSF  and  intensity 
model  must  be  estimated  from  the  data.  In  the  case  of  a  true  point  source,  as  modeled  in 
the  simulated  data  described  in  the  next  section,  the  object  is  treated  as  a  Kronecker  delta 
function,  S.  A  point  source  is  not  always  available  and  in  those  cases  the  object,  O,  must 
be  deconvolved  from  the  data  to  estimate  the  PSF.  In  addition,  the  data  contains  a  flat 
background  that  must  be  estimated. 

Estimates  of  the  intensity  model  and  PSF  are  made  from  the  data  using  an  EM  algorithm 
similar  to  one  developed  by  Schulz  [57].  The  difference  between  the  data  used  in  Schulz’s  blind 
deconvolution  and  the  data  used  for  the  following  derivation  is  that  the  data  includes  a  flat 
background.  The  intensity  model,  /,  as  a  function  of  pixel  coordinates  and  wavefront  error  is 

l(m,W)  =  0-h(m,W)  +  B ,  (5.5) 

where  the  estimate  of  the  background,  B, 

B  =  median  [<7  (m)  V/n  e[l,Md]],  (5.6) 

is  added  to  the  scaled  PSF.  Then  the  scale  factor,  6,  is  the  number  of  estimated  photons 
in  the  measured  point  source 

Mcl 

0  =  2^d(m)-B  ,  (5.7) 

m= 1 
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and  Md  is  the  number  of  pixel  coordinates  in  the  PSF  window  [18].  The  EM  algorithm 


used  to  estimate  the  PSF  begins  by  hypothesizing  the  data,  d  (m),  to  be 

d  (m)  =  dx  (m)  +  d2  (m)  (5.8) 


where  dx  (m)  and  d2(m ) 


are  independent  Poison  random  variables  with  mean  that  are  a 


functions  of  the  object  intensity,  0, 


dx{m)  =  9-h(m ), 


(5.9) 


and  of  the  background,  B, 


d2  (in) 


=  B. 


(5.10) 


Since  that  the  hypothesized  data  sets  are  assumed  to  be  independent  of  one  another  and 
Poison  distributed,  the  joint  PDF  is 

0-h(m)d'[m) eeh(m)  B‘h(m)e~B 


Md 

p(dvd2Vmel:Md)  =  Y}- 

x=\ 


dx(m)\  d2(m)\ 


and  the  log-likelihood  function  is 


Md 


L(0,h(m),B)  =  Yd,  (m)\n(d  ■  h  (m)) 

,n=\ 

-9  ■  h  (m)  -  d2  (m)ln  (5)  -  B. 


The  E-step  in  the  EM  algorithm  is  found  to  be 
0  =  £ 


M,, 


L[d,h(w)-,B)  |  d (m),9,h°ld  (w),5  -  A,Yji (m) 

m= 1 

9  ■  h  (m)~\-  9  ■  h  (m) 


tk 9-hold(m)-d(m ) 

=  2_j— - — In 


9-hold(m)  +  B 
B  -d  (m) 


9  ■  hold  +  B 


Md 

In  (5)-5-t£/7(/u), 

m= 1 


(5.11) 


(5.12) 


(5.13) 
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where  X  is  a  Lagrange  multiplier  [58].  The  Lagrange  multiplier  is  introduced  to  constrain 
the  estimation  of  the  PSF  to  sum  to  a  value  of  one.  Then  the  M-step  for  the  PSF  and 
Lagrange  multiplier  are  found  by  first  by  differentiating  0  with  respect  to  h  (m0) : 


50 

dh  (  mt) ) 


Md 

L 

m=l 


6-hoId  (m)-d\ 

9-hold  (m)  +  B 

h  (in) 

-e-A i. 


Therefore  the  PSF  update  for  each  pixel,  m0,  is  found  to  be 


(5.14) 


A”K)  = 


6  -h°,d  ( mQ )  ■  d  ( m0 ) 


d-hold(m0)  +  B\(0  +  A l) 


(5.15) 


and  the  Lagrange  multiplier  updates  to  be 


~d-hold(m) 


+  B 


(5.16) 


since 


Md 

'£Jh(m)  =  l.  (5.17) 

m= 1 


The  EM  algorithm  is  run  until  the  sum  of  squared  error  between  d  and  /  equals  the 
estimated  local  noise  variance  as  the  sum  of  d  and  at  that  point  the  PSF  is  estimated  to  be 

h(m)  =  h"(mll).  (5.18) 


5.2  Direct  Search  LS  vs.  E-Field  based  Estimation 

Zernike  coefficient  estimation  algorithms  that  rely  only  on  gradient  descent 
techniques  for  minimizing  the  sum  of  squared  error  (SSE)  between  the  modeled  and 
measured  intensity  patterns  have  the  problem  of  getting  trapped  in  local  minima  [35,  59]. 
A  new  Zemike  coefficient  technique  introduced  in  this  chapter  that  maximizes  the 
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correlation  between  an  estimated  electric  field  and  a  modeled  electric  field  also  has  the 


problem  of  getting  trapped  in  local  maxima.  However,  by  using  the  two  Zemike 
estimation  techniques  in  conjunction,  accurate  estimates  can  be  found  for  low  order 
Zernike  coefficients  (Z4-Z11)  that  correspond  to  the  aberrations  that  does  not  suffer  from 
local  minima/maxima. 

5.2.1  E-Field  versus  Intensity  Pattern  Correlation 

The  direct  search  LS  phase  retrieval  technique  uses  a  modeled  PSF  and  the 
measured  intensity  of  a  point  source  to  estimate  the  Zemike  coefficients.  However,  the 
modeled  electric  field  described  in  Eq.  (5.3)  can  also  be  used  to  estimate  Zemike 
coefficients  from  estimates  of  electric  field  in  the  detector  plane.  The  advantage  of  using 
the  electric  field  to  estimate  the  Zemike  coefficients  is  better  understood  by  quantifying 
and  comparing  examples  of  the  pair  wise  correlation  between  the  electric  field  patterns 
with  the  pair  wise  correlation  of  the  intensity  patterns  produced  by  the  Zernike 
polynomials  and  corresponding  Zernike  coefficient  shown  in  Figure  28.  In  this 
comparison,  the  modeled  electric  field  phase  is  completely  known  a  priori.  However,  in 
the  new  phase  retrieval  algorithm  presented  in  Section  5.2.4,  the  electric  field  is  estimated 
using  GS  and  may  not  produce  an  electric  field  that  is  as  accurate  as  if  was  directly 
measured  from  the  point  source.  So  this  demonstration  may  not  predict  how  well  the 
proposed  estimator  will  work,  but  instead  motivates  the  design  of  an  algorithm  that 
attempts  to  use  estimated  electric  fields  in  additional  to  intensity  patterns  for  Zernike 
coefficient  estimation. 
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Examples  of  the  numerical  analysis  results  displayed  as  surface  plots  are 
produced  using  the  Zernike  Polynomials  for  defocus,  (j)A,  and  spherical  aberrations,  (f)n. 

The  electric  field  associated  with  each  scaled  Zernike  polynomial  is  produced  using  Eqs. 
(5. l)-(5.3)  such  that 

H  (m,  Z.  .)  =  J7{a(u1)  exp  [j  ■  ZtJ  ■  (mj  )]},  (5.1 9) 

where  A(ux)  is  the  circular  aperture  function  shown  in  Figure  28. 


xIO'3 


-2  0  2 

Pupil  Plane  x-Position  [m]  x  ^q-3 


Figure  28.  Image  of  the  aperture  function,  A(ut ),  on  a  128  by  128  grid  for  the  experimental 
setup  described  in  this  chapter. 


Then  the  pair  wise  complex  correlation  coefficient  of  the  electric  field  pattern  in  the 
image  plane  is  computed  as 


RE 


P„(Z^Z,)  = 


(  mj  A/2  f  Md 

YH{m,Zi )  •  H  (m,  Z  )’  £//(«,  Z. )  •  H  (m,Z. ) 

U=i  y  U=i  y 


\/2 


(5.20) 
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where  RE  is  the  real  part  of  the  complex  value  and  the  pair  wise  nonnalized  correlation 
coefficient  of  the  intensity  pattern  in  the  image  plane  as 


ft(z,.z> 


^h(m,Z  ) ■  li(m,Z  j 

_ m=l _ 

Md  aX  (mz 

Tjh{m,Z,)2  •  Yjh{m,Zj) 


(5.21) 


The  correlation  coefficient  functions  measures  the  similarity  between  the  patterns 
on  a  scale  from  -1  to  1.  A  value  of  1  indicates  perfect  correlation,  a  value  of  0  indicates 
uncorrelated  and  a  value  of  -1  indicates  inverse  correlation  of  the  two  patterns  [56]. 

Figure  29  (a)  illustrates  the  difficulties  in  accurately  estimating  the  Zemike  coefficient  for 
defocus,  Z4,  in  an  optical  system  near  focus  with  only  the  intensity  data  due  to  the  high 
correlation  between  the  intensity  patterns  produced  by  different  amounts  of  defocus.  In 
contrast,  Figure  29  (b)  shows  that  the  E-field  pattern  produced  by  defocus  is  only 
perfectly  correlated  with  E-field  patterns  produced  by  one  with  equal  amount  of  defocus 
and  the  correlation  goes  down  dramatically  as  a  function  of  the  difference  in  defocus  that 
produces  the  two  E-field  patterns.  Again  in  Figure  29  (c)  the  intensity  patterns  produced 
by  spherical  error  and  defocus  are  highly  correlated,  where  as  the  electric  field  patterns 
are  significantly  less  so  as  shown  in  Figure  29  (d).  The  lower  correlation  between  the  E- 
field  patterns  makes  it  easier  to  distinguish  between  the  pattern  produced  by  defocus  from 
the  pattern  produced  by  spherical  error.  The  difference  in  correlation  in  the  E-field  space 
versus  the  intensity  space  is  considered  the  reason  why  the  electric  field  estimates  of  the 
Zernike  coefficients  differ  from  the  intensity  based  LS  phase  retrieval  estimates.  Those 


differences  enable  the  algorithm  described  later  in  this  paper  that  uses  both  E-field  and 
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intensity  estimates  of  Zernike  coefficients  iteratively  to  converge  on  the  correct  Zemike 
coefficients. 


Figure  29.  The  pair  wise  correlation  between  two  different  (a)  intensity  patterns,  ph,  and  (b) 
electric  field  patterns,  pH,  produced  in  the  image  plane  by  varying  the  Zemike  coefficients  for 
defocus,  Z4,  independently.  The  pair  wise  correlation  between  two  different  (c)  intensity  patterns 
and  (d)  electric  field  patterns  produced  in  the  image  plane  by  varying  the  Zemike  coefficient  for 
defocus  and  spherical  error,  Zn 

5.2.2  Least  Squares  Zernike  Coefficient  Estimation 

In  theory,  intensity  LS  phase  retrieval  could  be  accomplished  using  a  grid  search 
method  to  estimate  the  Zemike  coefficients,  but  because  of  the  number  of  parameters 
required  to  fonn  the  grid,  it  becomes  computationally  challenging  [38].  Therefore,  direct 
search  or  gradient  search  intensity  LS  methods  are  used  to  estimate  the  Zemike 
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coefficients  [60],  However,  these  search  methods  can  become  trapped  in  local  minima, 
but  the  minima  can  be  escaped  by  defocus ing  the  optical  system  [35]. 

The  intensity  based  LS  phase  retrieval  technique  uses  the  modeled  PSF,  h 

from  Eq.  (5.4)  and  measured  intensity  data,  d  ( m ) ,  of  a  point  source  to  estimate  the 

Zemike  coefficients.  The  PSF  is  estimated  from  the  data,  h  ( m ) ,  using  Eq.  (5.18)  and 

the  model  PSF  is  formed  as  a  function  of  wavefront  error  using  Eqs.  (5.1)-(5.4).  Then  the 
sum  of  squared  errors  (SSE),  Q,  between  the  modeled  and  the  measured  intensity  is 
detennined  by 

Md  ^  2 

Q  h{m),h{m,W^  =  ]T^/?(m)-/?(m,ffi)  .  (5.22) 

m= 1 

To  minimize  the  SSE,  the  intensity  LS  phase  retrieval  begins  with  an  initial  wavefront, 
W°,d ,  where 

WM  (ut  )  =  Z2-</,2(ui)  +  ...  +  Zn-<f,n(u1),  (5.23) 

which  for  the  application  of  focus  and  alignment  of  three  mirror  telescopes  can  be  limited 
to  only  Zemike  terms  2-11  [17].  The  wavefront  is  adjusted  by  adding  or  subtracting 
scaled  Zemike  polynomials  such  that  new  wavefronts  are  formed 

H'“ew(u,)=W°“(ul)  +  Az-<kM 

■  (5.24) 

C"(",)  =  r”“(“1)  +  az-4(",) 

W^(u,)  =  W"u(u,)-^  ■*,(«,) 

where  Az  is  the  scale  factor  that  determines  how  fine  of  an  increment  will  be  searched. 
Each  new  wavefront  error  results  in  the  following  new  modeled  PSFs 
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(5.25) 


h  [m,  )]  =  |  f{A(Ul ) exp [j  ■  W^(u, )]}|2 

h  [m,  W^fu, )]  =  |  J  [A{ux )  exp  [j  ■  W»m\u,  )]}|2 

h  [m,  )]  =  |  J  [A(u, )  exp  [j  ■  W^Xih)]}  f 

h  [m,  W^(u, )]  =  |  J{a[ux  )  exp  [j  ■  W^{u^\ . 

The  20  new  PSFs  in  Eqs.  (5.25)  are  each  substituted  into  Eq.  (5.22)  to  form  a  vector  of 
SSE  values 

~Q{hrn),h\m,W£m(uS§ 

Q{km),h\m,W^(uS§ 

Q=  \  ■  (5-26) 

Q  \h(m), h  [/w,  FT1JeM'(«j)]} 

Q{km),h[m,WiN;w(ux)^ 

The  wavefront  error  that  generates  the  minimum  SSE  value  in  Q  becomes  the  initial 
wavefront  error,  Wold{u1),  for  the  next  iteration  and  this  direct  search  process  continues 

until  the  wavefront  error  reaches  local  minima  [60].  Since  the  Zernike  coefficients  are 
used  to  fonn  each  new  wavefront  error  they  are  known  when  the  algorithm  ends. 

5.2.3  E-Field  Zernike  Coefficient  Estimation 

The  Zernike  estimation  method  using  the  estimated  electric  field  is  very  similar  to 
the  intensity  based  least  squares  method  described  in  the  last  section.  An  estimation  of 

the  PSF  from  the  data,  h(m ) ,  is  input  into  the  GS  algorithm  to  estimate  the  electric  field, 
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H[m ),  as  described  in  Section  2.2.2.  The  correlation  between  H (w)  and  the  modeled 
electric  field,  H(m,W),  is  found  as 


p(H(m),H{m,W )) 


RE 


Ma 


m—\ 


f  M,, 


X  1/ 
\/2 


j  (/n)  ■  H (; m )  j  (/ra, W) ' H (m> W) 


\m= 1  y  \w=l  y 

Again,  the  model  wavefronts  are  formed  by  adding  or  subtracting  scaled  Zemike 
polynomials  using  Eq.  (5 .24)  and  the  resulting  E-field  models  are 


(5.27) 


H  [m,  Wp(u, )]  =  J  {A  («, )  exp  [7  ■  (f2*>,)] 
H  \  m,  Wp(u,  )1  =  f(A(Ul  (expfy  •  Wp(ut) 


H  [m,  W“r(u, )]  =  ?{A  («,  )exp[y  •  («,)]} 

H  [m,  Wp(u, )]  =  f{A(ut  )exp[y  •  Wp(ut)]\. 


(5.28) 


The  20  new  electric  field  models  from  Eqs.  (5.28)  are  each  substituted  into  Eq.  (5.27)  to 
form  a  vector  of  correlation  values 


p{H{m),H[m,W^{uS§ 

p{H(m),H[m,W^(u^ 


(5.29) 


The  wavefront  error  that  corresponds  to  the  maximum  correlation  value  in  p 
becomes  the  initial  wavefront  error  for  the  next  iteration  and  this  direct  search  process 
continues  until  the  wavefront  error  reaches  a  local  maxima.  The  accuracy  of  the 
estimated  Zernike  coefficients  depends  on  how  the  algorithm  is  initialized  for  two 
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reasons.  The  first  reason  is  that  the  E-field  is  estimated  using  GS  and  not  a  measured  E- 
field.  The  second  reason  is  that  the  method  of  direct  search  is  used,  which  can  suffer 
from  getting  trapped  in  local  maxima.  The  next  section  discusses  an  iterative  process  that 
prevents  the  new  E-field  based  phase  retrieval  algorithm  from  estimating  the  wrong 
Zernike  coefficients  by  using  both  the  intensity  based  LS  and  E-field  correlation  direct 
search  methods  together. 

5.2.4  New  Phase  Retrieval  Algorithm 

The  new  phase  retrieval  algorithm  involves  an  iterative  process  that  includes  the 
E-field  estimation  from  the  data,  E-field  based  Zernike  coefficient  estimation,  and  direct 
search  LS  Zernike  coefficient  estimation  as  depicted  in  Figure  30.  The  PSF  is  estimated 
from  a  128  by  128  pixel  window  of  the  intensity  data  using  the  EM  algorithm  described 
in  Section  5.1.2  and  is  input  into  the  GS  algorithm  to  estimate  the  electric  field.  The 
estimated  electric  field  is  correlated  with  the  models  of  the  electric  field  to  estimate 
Zernike  coefficient  as  described  in  Section  5.2.3.  The  output  coefficients  become  the 
initial  conditions  for  the  intensity  based  LS  Zernike  coefficient  estimation  detailed  in 
Section  5.2.2.  If  the  SSE  value,  Q,  is  not  below  the  set  threshold,  a  new  wavefront  error 
is  modeled  using  Eq.  (5.1)  to  reinitialized  the  GS  algorithm.  The  process  is  repeated  until 
the  SSE  is  less  than  10'4  which  produces  accurate  estimates  of  the  Zernike  coefficients  in 
simulation  and  is  achieved. 
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Figure  30.  Block  diagram  new  electric  field  based  phase  retrieval  algorithm  for  estimating 
Zemike  coefficients. 


5.2.5  Phase  Retrieval  Simulation 

A  comparison  of  the  performance  of  the  new  phase  retrieval  algorithm  and  only 
the  intensity  based  LS  subroutine  were  simulated  to  illustrate  the  performance  advantages 
of  the  new  method.  The  following  intensity  model  is  used  to  simulate  data 

d  (m,W)  =  0  ■  h  ( m,W )  +  B  +  n(m),  (5.30) 

where  d(m,W)  is  Poisson  distributed,  0  is  the  number  of  photons  in  the  point  source,  B  is 
the  flat  background  noise,  and  n(m)  is  added  noise.  From  the  data,  estimates  of  the  point 
source  intensities  were  made  by  using  Eq.  (5.7)  and  background  was  found  with  Eq.(5.6). 
The  EM  algorithm  then  was  used  to  estimate  the  PSF. 

A  sample  of  simulated  inputs  and  outputs  of  both  the  new  electric  field  based 
phase  retrieval  method  and  the  intensity  based  LS  phase  retrieval  algorithm  are  shown  in 
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Figure  3 1 .  The  input  wavefront  error,  W,  shown  in  Figure  3 1  (a)  was  generated  setting 
all  the  Zernike  coefficients  Z4-Z11  equal  to  one.  The  simulated  data  from  W  is  shown  in 
Figure  31  (b)  and  was  generated  by  adding  a  100  photon  background,  a  104  photon 
source,  and  shot  noise  using  Eq.  (5.30).  An  image  of  the  phase  retrieved  model  of  the 
electric  field,  H(m,W)  shown  in  Figure  3 1  (d),  illustrates  the  how  the  pattern  differs  from 
that  of  the  phase  retrieved  PSF,  h(m,W),  shown  in  Figure  3 1  (e).  A  comparison  of  the 
phase  screens  of  wavefront  error  in  the  pupil  in  Figure  3 1  (a),  (c),  and  (f)  illustrates  the 
improved  accuracy  of  the  new  electric  field  based  phase  retrieval  algorithm  over  just 
using  the  intensity  based  direct  search  LS  algorithm.  Likewise,  by  comparing  the  PSF 
detennined  with  the  new  electric  field  based  phase  retrieval  algorithm,  Figure  3 1  (e),  and 
only  the  intensity  based  direct  search  LS  algorithm,  Figure  3 1  (g),  with  the  simulated  data 
in  Figure  3 1  (b)  it  is  evident  that  the  new  algorithm  is  performing  better  by  escaping  the 
local  minima/maxima. 
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Figure  31.  Phase  retrieval  simulation  examples  (a)  input  wavefront  in  pupil  (b)  simulated  data  in 
image  plane  (c)  estimated  wavefront  error  in  pupil  using  new  electric  field  based  phase  retrieval 
algorithm  (d)  real  part  of  the  electric  field  pattern  in  image  plane  determined  using  new  electric 
field  based  phase  retrieval  algorithm  (e)  PSF  in  image  plane  determined  using  new  electric  field 
based  phase  retrieval  algorithm  (f)  estimated  wavefront  error  in  pupil  using  only  the  intensity 
based  direct  search  LS  algorithm  (g)  PSF  in  image  plane  determined  using  only  the  intensity 
based  direct  search  LS  algorithm. 


To  evaluate  the  performance  of  the  phase  retrieval  algorithm,  100  combinations  of 
random  low  order  Zemike  coefficients  were  used  to  produce  the  simulated  wavefront 
error  where 

[Z2  Z3  ...  Zn]  =  unif(- 1,1),  (5.31) 

where  the  unif  function  is  the  uniform  distribution.  Each  of  the  randomly  generated 
wavefront  errors  has  an  associated  an  intensity  model  that  has  noise  added  to  produce  100 
trials,  T,  of  data  using  Eq.  (5.30).  The  Zernike  coefficients  are  estimated  using  both  the 
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intensity  based  direct  search  LS  algorithm  and  the  new  electric  field  based  phase  retrieval 
algorithm  starting  with  the  Zernike  coefficients  initialized  to  zero.  The  results  of  each 
algorithm’s  Zernike  coefficient  estimates,  Zt,  are  plotted  in  Figure  32  and  are  in  terms  of 

the  difference,  A,  between  the  absolute  value  input  of  the  known  Zernike,  Z(,  and  the 
average  estimated  Zernike  coefficients 


A=  Z.J 


Z*. 
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where  i  e  [4,...,1  lj . 


(5.32) 


The  mean  A  is  computed  for  all  the  random  combinations,  C,  of  Zernike  coefficients  used 
to  fonn  the  simulated  wavefront  error 
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ZX 


A  --£=! 

Mean 
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and  the  standard  deviation  of  A  is 


100 


,C=1 


A  stD  =  \SL  {^C-^Mean) 


/100 


(5.33) 


(5.34) 
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□  Intensity  Based 
x  E-Field  Based 


x  E-Field  Based 


Figure  32.  Difference  between  the  Zemike  coefficients,  Z,,  simulated  and  phased  retrieved  using 
(a)  only  the  intensity  based  direct  search  LS  algorithm  versus  new  electric  field  based  phase 
retrieval  algorithm  and  (b)  expanded  plot  of  new  electric  field  based  phase  retrieval  algorithm 
results.  The  error  bars  corresponded  to  the  standard  deviation  of  A  are  computed  via  Eq  (5.34). 


Using  only  the  intensity  based  direct  search  LS  algorithm  produced  the  biased 
estimates  of  the  Zernike  coefficients  due  to  the  problem  of  getting  trapped  in  local 
minima  as  shown  in  Figure  32  (a).  In  contrast,  the  estimates  of  the  Zemike  coefficients 
using  the  new  electric  field  based  phase  retrieval  algorithm  expanded  on  the  plot  in 
Figure  32  (b)  are  relatively  unbiased  because  the  local  minima  are  avoided  due  to  the 
iterative  process  outlined  in  Section  5.2.4. 


5.3  Laboratory  Demonstrations 

To  investigate  the  performance  of  the  phase  retrieval  algorithm  on  measured  data 
two  different  controlled  aberrations  were  generated  in  the  laboratory.  The  first  was 
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defocus,  which  is  simple  to  generate,  control  and  quantify  by  moving  the  image  plane 
with  respect  to  the  detector.  The  second  was  astigmatism  generated  by  tilting  the  lens 
with  respect  to  the  optic  axis  which  demonstrates  perfonnance  on  another  aberration  but 
is  more  difficult  to  control. 


5.3.1  Defocus 

The  defocus  demonstration  used  the  setup  depicted  in  Figure  33.  The  object,  O, 
was  fonned  by  a  LED  with  a  wavelength,  X,  of  648  mn  illuminating  a  200  pm  pinhole. 
The  object  is  imaged  by  a  single  lens  with  a  focal  length,  /  ,  of  17.5  cm.  The  focused 


image  was  fonned  in  image  space  according  to  the  thin  lens  equation  [5] 

i-_L  _L 
7 =^+V 


(5.35) 


where  SQ  is  the  distance  from  pin  hole  to  the  aperture  and  S,  is  the  image  location.  The 


aperture  is  defined  by  an  aperture  stop  with  diameter,  D,  of  3.62  mm  to  ensure  that  the 
Nyquist  sampling  was  achieved  with  the  16  pm  pixel  pitch  of  the  camera’s  detector. 
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Figure  33.  Defocus  demonstration  setup. 


Defocus  was  generated  by  shifting  the  lens  mounted  on  a  translation  stage  in 
order  to  defocus  the  image  by  creating  a  distance  between  the  image  location  and  the 
detector  distance,  SD,  from  the  lens 

A  focus  =  S;  ~  SD’  (5-36) 

according  to  Table  4. 


Table  4.  Defocus  Demonstration  Parameters 


So  (cm) 

Si  (cm) 

^Focus  (cm) 

188.2 

19.8 

-0.5 

189.2 

19.3 

0 

190.2 

18.8 

0.5 

191.2 

18.3 

1 

192.2 

17.8 

1.5 
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The  theoretical  wavefront  error  that  was  generated  from  these  shifts  is  computed  and 
decomposed  into  the  Zernike  coefficients,  Zi-Zn  [31]. 

Figure  34  summarizes  the  results  of  the  defocus  demonstration.  The  first  column 
is  the  measured  shift  in  the  lens,  which  is  equal  to  the  focus.  The  second  column  is  that 
intensity  data  that  was  recorded  by  the  camera  corresponding  to  the  lens  position.  As 
expected  with  more  defocus,  the  data  becomes  spread  out.  The  PSF  column  shows  the 
data  after  it  is  blind  deconcolved  from  the  pin  hole  object  to  produce  the  PSF,  since  the 
pinhole  was  not  small  enough  to  represent  a  true  point  source  [16].  The  blind 
deconvolution  is  achieved  through  a  two-step  process.  First  the  blurred  image  of  the 
pinhole  is  estimated  using  the  EM  algorithm  presented  in  Section  5.1.2.  This  has  the 
effect  of  producing  a  background  removed  image  of  the  pinhole  convolved  with  the  point 
spread  function.  Then  Schulz’s  blind  deconvolution  algorithm  (BDA)  estimated  the 
pinhole  shape  and  the  point  spread  function.  A  circular  function  is  used  to  initialize  the 
image  estimate  and  phase  function  and  two  waves  of  defocus  are  used  to  initialize  the 
filed  in  the  pupil.  The  pupil  transmittance  function  is  circular  with  a  diameter  of  3.62 
mm  [57].  An  initial  estimate  of  the  pinhole  size  is  input  into  the  BDA,  in  terms  of  the 
number  of  pixels  in  object  space,  along  with  the  nonnalized  data  output  from  the  EM 
algorithm.  The  BDA  is  run  for  100  iterations  in  which  it  converges  on  a  stable  estimate 
of  the  size  and  shape  of  the  pinhole  while  simultaneously  estimating  the  PSF.  Phase 
retrieval  is  conducted  on  the  PSF  to  estimate  Zemike  coefficients  (Z2-Z11)  and  produce 
the  modeled  PSF.  The  estimated  coefficient  for  defocus  is  listed  next  to  the  theoretical 
value  that  is  computed  using  the  wave  optics  calculations  discussed  in  the  appendix.  The 
small  difference  between  the  theoretical  value  for  Z4  and  the  estimate  are  within  the 
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combined  error  of  our  lens  translation  stage  accuracy  of  ~1  mm  and  the  estimate  error 
from  noise.  One  important  point  to  note  is  that  only  defocus  is  experimentally  generated, 
but  the  phase  retrieval  algorithm  attempted  to  estimate  10  different  Zemike  coefficients, 
yet  it  still  correctly  attributed  the  wavefront  error  to  the  aberration  of  defocus. 


Lens 

Translation 

A  =  -.5cm 


A  =  0  cm 


Data  PSF  Model 


A  =  0.5  cm 


A  =  1  cm 


A  =  1.5  cm 


|Z4|  Theory 
±.13  [waves] 

|Z4|  Est 
±.03[waves] 

.60 

.75 

.00 

.00 

.61 

.60 

1.26 

1.15 

1.95 

2.10 

Figure  34.  Results  of  defocus  demonstration. 


5.3.2  Astigmatism 

The  Astigmatism  demonstration  uses  the  set  up  depicted  in  Figure  35.  The  object, 
O,  is  formed  with  the  same  pinhole  and  LED  as  the  defocus  demonstration  by  a  lens  with 
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a  focal  length,  f  of  50.0  cm.  The  camera  detector  is  placed  at  Si  and  the  lens  is  tilted  with 
respect  to  the  optic  axis  in  order  to  generate  astigmatism. 


Figure  36  shows  the  results  of  the  astigmatism  demonstration.  The  first  column  is 
the  camera  recorded  data  taken  in  and  out  of  focus.  The  data  is  processed  with  the  blind 
deconvolution  algorithm  to  as  with  the  defocus  experiment  produce  the  PSF  column 
Figure  36  (b),  (e),  and  (h).  Then  the  Zemike  model  is  phase  retrieved  from  the  focused 
PSF  in  Figure  36  (b).  The  Zemike  model  for  the  defocus  spots  in  Figure  36  (f)  and  (i)  are 
fonned  by  using  the  model  terms  Z5-Z11  from  the  phase  retrieval  of  the  focus  spot  and 
then  fit  using  Z2-Z4  as  free  parameters.  Comparison  of  the  PSF  in  Figure  36  (e)  and  (h) 
to  the  modeled  PSF  in  Figure  36  (f)  and  (i)  illustrates  that  an  accurate  amount  of 
astigmatism  in  the  optical  system  is  phase  retrieved  from  the  focused  spot. 
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Table  5  lists  the  Zernike  coefficients  (Z5-Z11)  that  are  retrieved  from  the  focused 
PSF  shown  in  Figure  36  (b).  The  coefficients  are  used  to  produce  the  model  displayed  as 
color  maps  in  Figure  36  (c),  (f)  and  (i).  Note  that  the  astigmatism  term,  Z6,  quantifies  the 
estimated  amount  of  astigmatism  generated  in  the  tilted  lens  demonstration.  The  in¬ 
focus,  negative-focus,  and  positive-focus  columns  correspond  to  the  position  of  the  CCD 
relative  to  focus  and  the  resulting  value  of  the  phase  retrieved  Zernike  coefficient  for 
defocus.  The  theoretical  column  is  computed  using  the  wave  optics  calculations 
discussed  in  the  appendix  and  tilting  the  lens  at  an  angle  of  7.2  deg.  As  predicted  by 
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theory,  tilting  the  lens  generates  a  measurable  amount  of  astigmatism  and  the  new  phase 
retrieval  method  is  able  to  quantify  that  aberration  while  jointly  estimating  the  other 
parameters  at  focus.  The  small  biases  of  some  of  the  estimated  Zemike  coefficients  from 
the  theoretical  values  (<0.15  waves)  listed  in  Table  5  are  likely  due  to  noise  in  the  data. 


Table  5.  Astigmatism  Demonstration  Zemike  Coefficient  Estimates 


Zemike 

Coefficient 

In- 

Focus 

Negative- 

Focus 

Positive- 

Focus 

In-Focus 

Theory 

z4 

0 

-4.6 

3.9 

0 

z5 

0.15 

0.15 

0.15 

0 

Z6 

-0.8 

-0.8 

-0.8 

-0.8 

z  7 

-0.05 

-0.05 

-0.05 

0.1 

Z8 

-0.1 

-0.1 

-0.1 

0 

Z9 

0.05 

0.05 

0.05 

0 

Zio 

0 

0 

0 

0 

Zn 

-0.1 

-0.1 

-0.1 

0 

5.4  Conclusions 

The  iterative  use  of  three  different  phase  retrieval  techniques:  Gerchberg-Saxton, 
estimated  electric  field  correlation,  and  intensity  LS  produces  accurate  estimates  of  the 
magnitude  of  low  order  Zemike  coefficients  of  an  optical  system  from  focused  point 
source  images.  Each  of  these  algorithms  can  estimate  the  Zernike  coefficients  alone,  but 
have  demonstrated  biases  that  are  related  to  how  they  are  initialized.  By  using  the  three 
methods  together,  the  need  to  guess  the  starting  wavefront  error  has  not  been  necessary 
for  the  algorithm  to  accurately  estimate  the  Zemike  coefficients.  In  addition,  with  this 
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approach  the  limitations  posed  by  the  wrapped  2-D  phase  produced  by  the  GS  algorithm 
are  avoided.  Differences  in  the  electric  field  pattern  and  intensity  pattern  enable  the 
algorithm  to  escape  from  local  minima/maxima  and  converge  on  the  correct  solution. 
Simulations  and  laboratory  results  demonstrate  the  perfonnance  of  the  algorithm. 

In  the  future,  the  new  iterative  algorithm  could  potentially  be  improved.  One 
concept  is  to  use  other  optical  constraints  on  the  GS  portion  of  the  algorithm  (like  in  the 
hybrid  input  output  algorithm)  to  improve  Zernike  estimates  and  decrease  convergence 
time  [61].  Another  idea  is  to  use  a  gradient  search  method  (such  as  the  Levine  Marquette 
method)  to  speed  up  the  LS  portion  of  the  algorithm. 

The  ability  to  accurately  estimate  the  Zemike  coefficients  while  the  optics  are 
near  focus  should  be  useful  in  the  alignment  of  large  three  mirror  telescopes  like  the  SST 
and  LSST  because  the  aberration  will  not  be  reintroduced  from  movements  of  the 
secondary  and  tertiary  mirrors  required  for  focusing.  Using  this  phase  retrieval  method 
for  JWST  would  simplify  their  optical  design  and  reduce  engineering  risk  by  removing 
the  need  for  different  focal  length  lenses  to  generate  different  defocus  values  for  focus- 
diverse  phase  retrieval  [45,  47,  48].  In  addition,  the  new  phase  retrieval  algorithm 
proposed  in  this  chapter  has  the  advantage  over  NASA’s  hybrid  diversity  algorithm 
(HD A)  planned  for  use  on  JWST  because  the  HDA  uses  not  only  phase  diversity  but  also 
has  to  unwrap  the  pupil  phase  as  its  estimates  are  based  on  the  phase  output  from  GS 
[47]. 
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VI.  Conclusions 


Substantial  efforts  have  been  made  by  astronomers  to  maximize  the  ability  to 
detect  space  objects  with  optical  telescopes,  but  there  is  a  continued  need  to  optimize  the 
process  for  both  military  and  scientific  applications.  Improved  image  processing  affords 
a  way  to  improve  space  object  detection  without  the  need  to  invest  in  more  costly 
hardware  like  space  based  telescopes  or  larger  aperture  ground  based  telescopes.  An 
evaluation  of  current  image  processing  techniques  was  conducted  to  identify  how  the 
SST’s  image  data  could  be  better  processed.  The  two  areas  uncovered  were  phase 
retrieval  methods  and  detection  algorithms,  which  are  summarized  respectively  in 
Section  2.2  and  Section  2.4. 

A  thorough  literature  review  of  existing  phase  retrieval  techniques  and  detectors 
available  for  use  on  the  SST  are  summarized  in  Section  2.5.  The  primary  limitation  of 
present  phase  retrieval  techniques  is  that  they  require  phase  diversity  that  is  typically 
induced  by  defocusing  the  telescope  as  discussed  in  Section  2.2.  Defocusing  the 
telescope  is  undesirable  because  of  complicated  defocusing  and  alignment  procedures  of 
three  mirror  telescopes.  In  terms  of  detection  improvement,  one  possible  detection 
algorithm  identified  in  literature  that  could  be  implemented  on  the  SST  in  place  of  its 
baseline  point  detector  is  the  correlator  discussed  in  Section  2.4.  However,  previous 
AFIT  research  has  shown  that  the  correlator  suffers  from  aliasing  when  the  data  is 
undersampled  [43]. 

The  image  processing  techniques  that  were  developed  in  the  pursuit  of  this 
research  have  demonstrated  the  potential  to  improve  the  SST  in  tenns  of  its  system 
perfonnance  metrics  described  in  Section  1 . 1  including:  metric  accuracy,  photometric 
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accuracy,  and  sensitivity.  The  new  phase  retrieval  techniques  that  produce  low  bias 
estimates  of  telescope  aberrations  in  tenns  of  Zemike  polynomials  are  illustrated  in 
Figures  18  and  32.  Based  on  the  results  shown,  they  could  improve  the  focus  and 
alignment  of  the  SST,  which  would  provide  immediate  improvement  in  the  telescope’s 
sensitivity  by  reducing  the  FWHM  of  its  PSF.  Even  more  noteworthy,  the  MHT  has 
superior  properties  that  could  be  used  as  the  kernel  of  an  entirely  new  image  processing 
scheme  that  would  significantly  outperfonn  the  baseline  algorithm.  One  of  those 
properties  is  sub-pixel  position  information  that  should  improve  metric  accuracy  as 
described  in  Section  4.6.  Another  property  of  the  MHT  is  the  linear  relation  between  its 
SNR  output  and  estimates  of  the  object’s  brightness,  which  has  demonstrated  improved 
photometric  accuracy  over  the  baseline  detector  as  shown  in  Figure  27.  The  most 
significant  improvement  in  telescope  perfonnance  accomplished  with  the  MHT  over  the 
exiting  software  comes  from  detection  sensitivity  gains  as  illustrated  in  the  GEO  eclipse 
experimental  analysis  results  plotted  in  Figure  26. 

6.1  Long  Exposure  Phase  Retrieval  Improvements 

Increased  understanding  of  the  long  exposure  phase  retrieval  problem  and  better 
ways  to  perform  it  have  been  developed  in  this  work.  One  fresh  insight  is  that  the 
presence  of  a  long  exposure  atmosphere  affects  the  CRLB  for  Zemike  coefficients  and 
that  the  bound  helps  predict  phase  retrieval  performance.  As  shown  in  Figure  12,  the 
long  exposure  atmosphere  increases  the  CRLB  for  Zernike  coefficient  for  defocus. 
Additionally,  Figure  16  shows  that  the  bound  predicts  the  general  trend  of  phase  retrieval 
performance  for  estimates  of  the  defocus  parameter.  Another  key  phase  retrieval  finding 
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is  that  unbiased  estimates  of  Zernike  coefficients  in  the  presence  of  an  atmosphere  can  be 
made  without  requiring  a  defocused  spot  to  generate  phase  diversity.  Chapter  III 
describes  a  grid  search  method  that  can  produce  unbiased  estimates  of  spherical  and 
defocus  coefficients  in  the  presence  of  a  long  exposure  atmosphere  as  seen  in  Figure  18. 

In  the  long  exposure  atmospheric  case,  the  limitations  of  differentiating  between 
the  contribution  of  telescope  aberrations  and  atmospheric  blurring  to  the  intensity  pattern 
are  overcome  by  jointly  estimating  the  atmospheric  seeing  parameter  and  Zernike 
coefficients  using  Eq.  (3.38).  The  method  works  because  it  uses  a  grid  search  method  to 
select  from  a  library  of  possible  PSFs  to  find  the  PSF  that  best  matches  the  measured 
intensity  pattern.  The  selected  PSF  corresponds  to  specific  Zernike  coefficients  and 
seeing  parameter  values.  This  method  is  useful  when  there  are  a  small  number  of 
coefficients  that  need  to  be  estimated.  This  is  the  case  for  the  SST  because  Z5-Zn  were 
previously  estimated  using  a  phase  diversity  techniques  and  only  ro  and  Z4  need  to  be 
estimated  near  focus  in  order  to  build  the  PSF  model  used  in  the  MHT  (see  Section  4.4). 

6.2  Telescope  Detection  Improvements 

The  implementation  of  a  MHT  by  the  SST  will  increase  its  probability  of 
detecting  space  objects  over  the  baseline  detector  through  three  improvements:  using 
multiple  pixels  to  compute  SNR,  mitigating  aliasing  effects,  and  computing  more 
accurate  window  statistics.  The  MHT  mitigates  the  combined  effects  of  atmospheric 
blurring  and  telescope  aberrations  that  cause  point  source  light  to  spread  across  multiple 
pixels.  This  is  done  by  including  the  PSF  in  the  MHT  through  Eq.  (4.13),  which 
computes  higher  SNR  values  than  the  baseline  detector  because  the  point  detector  only 
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uses  one  pixel  to  compute  SNR  through  Eq.  (2.32).  In  addition,  since  the  SST’s  data  is 
undersampled,  the  measured  PSF  changes  shape  when  the  light  is  not  centered  on  a  pixel 
and  the  MHT  is  able  to  mitigate  the  aliasing  effects  that  degrade  the  performance  of  both 
the  correlator  and  the  baseline  detector.  This  is  accomplished  by  shifting  position  on  the 
Nyquist  grid  as  discussed  in  Section  4.4.  In  addition,  the  computation  of  the  noise 
statistics  in  the  window  surrounding  the  pixel  being  tested  is  skewed  by  stars  crossing, 
cosmic  radiation  and  bad  pixels.  However,  by  removing  outlier  pixel  data,  more  accurate 
standard  deviations  are  computed  and  used  in  the  MHT  as  discussed  in  Section  4.5.  By 
combining  these  improvements  into  the  MHT  through  Eq.  (4.13)  and  Eq.  (4.18),  the 
MHT  demonstrates  improvement  in  Pd  of  over  50%  over  the  baseline  detector  based  on 
the  data  collected  in  the  ANIK-F1  GEO  eclipse  experiment  as  shown  in  Figure  26. 

6.3  Short  Exposure  Phase  Retrieval  Improvements 

In  future  the  SST  sites,  Zemike  coefficients  may  not  be  estimable  with  the 
currently  used  phase  diversity  techniques  due  to  poor  atmospheric  seeing  condition  or 
because  defocusing  the  telescope  may  be  undesirable.  Unfortunately,  the  long  exposure 
phase  retrieval  method  is  too  computationally  burdensome  for  estimating  large  numbers 
of  Zemike  coefficients.  However,  if  the  SST  upgrades  to  a  frame  transfer  camera,  the 
phase  retrieval  algorithm  described  in  chapter  V  can  produce  unbiased  estimates  of 
Zemike  coefficients  (Z4-Z11)  in  the  presence  of  the  short  exposure  atmosphere  as  shown 
in  Figure  32. 

The  primary  improvement  of  the  short  exposure  phase  retrieval  for  Zemike 
coefficients  over  existing  methods  is  that  it  does  not  require  defocusing  of  the  telescope 
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in  order  to  estimate  coefficients  thus  making  it  a  useful  new  tool  for  optical  diagnostics. 
The  short  exposure  phase  retrieval  of  Zernike  coefficients  from  a  focused  telescope  point 
source  intensity  pattern  using  only  an  intensity  based  least  squares  direct  search  approach 
is  challenging  because  of  the  correlation  between  intensity  patterns  produced  by  different 
aberrations  as  illustrated  in  Figures  29  (a)  and  (c).  In  contrast,  the  electric  field  patterns 
produced  by  the  same  aberrations  are  much  less  correlated  as  shown  in  Figures  29  (b)  and 
(d).  Therefore,  by  using  an  estimated  electric  field  in  conjunction  with  the  intensity 
pattern,  the  new  short  exposure  phase  retrieval  algorithm  is  able  to  produced  nearly 
unbiased  estimates  of  the  Zemike  coefficients  from  both  simulated  and  experimental  data 
as  discussed  in  Sections  5. 2-5. 3. 

6.4  Future  Work 

While  significant  strides  have  been  made  to  enhance  the  perfonnance  of  the  SST 
and  other  astronomical  telescopes,  many  other  interesting  and  relevant  questions  that 
could  be  answered  in  future  research  have  been  uncovered  including:  Is  it  better  to  use 
the  long  or  short  exposure  atmospheric  model  for  phase  retrieval  of  Zernike  coefficients? 
How  does  the  MHT  perform  as  a  function  of  atmospheric  seeing?  How  often  should 
phase  retrieval  be  conducted  to  mitigate  seeing  effects  on  the  detection  performance  of 
the  MHT?  How  much  degradation  in  MHT  performance  results  from  the  changes  in  PSF 
shape?  What  is  the  range  of  PSF  shapes  that  can  be  tolerated?  How  much  does  the 
performance  of  the  correlation  based  BHT  depend  on  the  PSF  shape?  How  is  the 
performance  of  the  MHT  affected  by  variations  in  the  PSF  across  the  FOV?  Will 
multiple  PSF  shapes  be  used  for  different  positions  on  the  field?  If  not,  what  is  the 
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degradation  in  perfonnance  if  a  single  PSF  is  used  for  the  field?  By  addressing  these 
questions  even  better  detection  performance  may  be  achieve  with  SST. 

6.5  Final  Observations 

Investigating  how  to  improve  image  processing  for  space  object  detection  has 
been  a  fruitful  area  of  research  and  a  productive  means  of  enhancing  ground  based 
telescope  performance.  While  better  ways  for  phase  retrieval  and  detection  have  been 
identified  in  this  dissertation,  it  has  also  opened  up  an  entire  new  area  for  future  research 
as  outlined  with  the  questions  in  the  previous  section.  Continued  pursuit  to  answer  these 
questions  and  the  new  ones  that  are  certain  to  arise  should  be  productive  both  from 
academic  and  DoD  points  of  view.  This  is  because  further  enhancement  to  image 
processing  of  the  SST  and  other  telescopes  contributing  to  the  space  surveillance  network 
has  the  potential  to  enhance  the  United  States’  SSA  capabilities. 
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VII.  Appendices 


The  following  three  sections  provide  addition  supporting  evidence  for  the  validity 
of  the  phase  retrieval  performance  described  in  Chapters  III  and  IV.  The  first  two 
sections  provide  empirical  evidence  that  the  SST  star  data  used  for  phase  retrieval  is 
modeled  with  an  appropriate  cumulative  distribution  function  (CDF)  and  that  the  pixels 
can  be  considered  statistically  independent.  The  last  section  covers  the  method  by  which 
the  theoretical  Zemike  coefficients  for  the  aberrations  produced  in  the  lens  experiment 
described  in  Chapter  IV  were  computed. 

A.l  The  SST’s  CCD  Noise  Statistics 

Raw  measurements  of  digital  counts  were  recorded  from  a  star  observed  with  the 
SST  over  multiple  frames.  An  empirical  quantile-quantile  (Q-Q)  plot  of  the  data  from  the 
pixel  the  star  was  centered  on  versus  a  Poisson  probability  mass  function  (PMF)  with  the 
same  mean  as  the  star  intensity  is  shown  in  Figure  37  (a)  [62].  By  inspecting  the  plot  it 
appears  that  the  Poisson  PMF  is  not  an  ideal  match  for  the  data  because  an  excellent  fit 
would  have  the  blue  data  marker  on  the  red  line  formed  by  the  ideal  distribution.  One 
possibility  for  the  disparity  between  the  data  distribution  and  the  Poisson  PMF  is  that  the 
data  is  quantized. 

Quantization  noise  occurs  in  the  analog-to-digital  conversion  process  from 
photons  received  to  digital  counts  if  it  is  not  a  one-to-one  conversion.  The  gain,  G,  that 
converts  between  the  total  number  photons  received  in  a  pixel  and  expected  digital  count 
is  represented  by 
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E[d(m)]  =  G  'y  OJijm  —x)  +  B 


(7.1) 


where,  E^d  (w)J  is  the  mean  number  of  digital  counts  in  the  pixel  of  interest,  m;  <9,  is  the 

intensity  of  the  space  object,  x  is  the  pixel  coordinates  in  the  window  of  interest; 
h(m  —  x)  is  the  PSF  and  B  is  the  background.  The  mean  number  of  photons  incident  on 
a  pixel 

£[cT(m)]  =  -x)  +  B.  (7.2) 

Since  gain  is  a  constant,  the  relation  between  the  digital  counts  and  the  photon  incented 
on  a  pixel,  d'(m),  is 

d  (m)  =  Gd'(m).  (7.3) 

Shot  noise  is  caused  by  the  fact  the  independent  photon  arrival  times  on  the  CCD 
for  each  pixel  in  the  image  can  be  considered  a  collection  of  Poisson  random  variables 
(RV),  therefore  d'(m)  has  a  Poisson  PMF  [2].  If  quantization  noise  and  shot  noise  are 

the  dominate  noises  in  d  (m),  then  d(m)/G  will  also  have  a  closer  match  to  Poisson 
statistics. 

To  characterize  the  CCD  and  determine  G  from  the  imagery  data,  the  following 
derivation  is  used.  For  a  given  pixel  the  mean  number  of  digital  counts  is 

d  =  fs[<i(/w)]  =  (7.4) 

and  the  second  moment  is 
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( d(m)-d )  =  £  d2  (m)  -2d(m)d  +  d2 


=  E[d2(m)]-2d2+d 2 
-  E^d2  (m)J- d2  (m). 


In  addition, 


d2(m)  =E  ( Gd'(m )) 


G2£ 


d'2  (m) 


(7.5) 


(7.6) 


Substituting  Eq.  (7.6)  and  Eq.  (7.4)  into  Eq.  (7.5) 

£[</2  (»;)]  -  d2  =  G2E[d'2  (m)]  -  G2E2  [d'(m)\ 

=  G2  {j E[d '2  (m)]  -  E2  [tf'O")]}  t7-1^ 

=  G2E[d'{m)\. 


because  for  a  Poisson  random  variable  (RV)  the  mean  is  equal  to  the  variance. 
Therefore,  G  can  be  estimated  from  the  data  by  dividing  Eq.  (7.7)  by  Eq.  (7.4) 

E[d2(m)\- d2  _G2E[d'(m)]_ 

E[d(m)]  _  GE[d'(m)\  ~ 


(7.8) 


Using  star  data  from  the  SST  the  gain  was  estimated  to  be  G  =  3. 

The  gain  estimate  was  used  to  remove  the  quantization  effect  from  the  raw  data 
shown  in  Figure  37  (a)  to  produce  the  adjusted  Q-Q  plot  in  Figure  37  (b).  The  plot  shows 
that  the  data  with  quantization  noise  removed  is  better  modeled  with  the  Poisson  PMF, 
because  the  magnitude  of  the  offsets  of  the  blue  data  points  from  the  ideal  distribution 
indicated  by  the  red  line  are  significantly  reduced.  This  analysis  provides  a  measure  of 
confidence  that  the  Poisson  PMF  should  be  used  for  the  derivation  of  the  LRT  detectors 
and  CRLB  for  estimates  of  Zernike  polynomial  coefficients. 


121 


Figure  37.  Q-Q  Plot  of  star  data  from  the  SST  images  versus  a  Poisson  distribution  (a)  raw  data 
(b)  Corrected  for  quantization 


A.2  The  SST’s  CCD  Pixel  Independence 

The  same  star  data  used  to  characterize  the  SST’s  CCD  noise  statistics  is  used  to 
evaluate  the  independence  of  the  data  in  each  pixel.  The  Spearman  coefficient  of  rank 
correlation,  p,  was  used  find  the  correlation  between  the  pixel  that  the  star  was  centered 
on  and  the  adjacent  pixels  as  shown  in  Figure  38  [62].  For  the  pixels  evaluated,  p  was 
much  closer  to  0  than  to  1  or  -1  indicating  that  the  pixels  are  independent.  Furthermore, 
the  p-value  is  larger  than  0.1  in  all  four  cases  indicating  that  the  correlation  is 
significantly  different  from  1 . 
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Figure  38.  Correlation  between  a  star’s  center  pixel  and  adjacent  pixels 


A.3  Aberration  Calculations  for  a  Single  Lens 

The  aberrations  produced  by  a  single  lens  are  calculated  using  a  wave  optics 
approach  to  verify  the  accuracy  of  the  phase  retrieval  results  in  the  laboratory 
demonstrations  presented  in  this  paper.  The  point  source  is  assumed  to  generate  a 
spherical  wavefront  that  is  propagated  to  the  lens.  The  lens  introduces  a  quadratic  phase 
transformation  and  then  the  spherical  wavefront  is  propagated  to  the  image  plane  [31]. 
The  combined  wavefront  error  from  the  two  propagations  and  the  lens  transformation  are 
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summed  up  in  the  pupil  and  then  decomposed  into  scaled  Zernike  polynomials  to 
detennine  the  theoretical  Zemike  coefficients. 

The  wave  optics  simulation  input  values  are  the  object  and  image  plane  distances 
from  the  lens,  S,  and  SQ  ,  the  focal  length  of  the  lens,  f,  the  wavelength  of  light,  X,  and 
the  diameter  of  the  entrance  pupil,  D.  Then  the  magnification  by  the  lens,  M,  is 
computed  as  [30] 


M  =  - 


S, 


(7.9) 


Two  variables  are  used,  one  is  the  amount  of  displacement  the  object  from  the  optic  axis, 
Ay,  and  the  other  the  is  distance  the  camera  is  from  image  plane,  A z. 

The  wave  optics  simulation  beings  with  the  generation  of  a  1024  by  1024  phase  screen, 
such  that  the  pixel  size  in  the  pupil,  dx,  is 


dx 


D 

1024' 


(7.10) 


Two  1024  by  1024  matrices  of  the  x  and  y  coordinates  are  generated,  x  and  V  ,  and  are 
set  such  that  the  center  pixel  in  the  grid  has  coordinate  (0,0)  and  have  the  pupil  spacing  of 
dx.  With  the  phase  screens  set  up,  the  object  is  modeled  as  a  spherical  wave,  Wv 

propagated  from  the  object  plane  to  the  lens  using  the  following  distance  fonnula  [3] 

W,  =2-x-f(v-Av)2  +S02Y2  /*■  (7-11) 


To  compute  the  lens  transfonnation,  the  radius  of  the  pupil,  r,  is  computed  as  a  function 
of  pixel  coordinates 


l  2  l)^2 

r  =  u +y_  ’ 


(7.12) 
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and  the  quadratic  phase  factor  for  the  lens,  A,  is  computed  as 


A  = 


-K  •  r 


(7.13) 


Next  the  wavefront  error  caused  by  propagation  from  the  lens  to  image  location  in  image 
plane,  W2,  is  computed  using  the  distance  fonnula 


W,  =2-tt- 


[y-M-Ay)  +(S'/-Az)2j  jx. 


(7.14) 


The  total  wavefront  in  the  pupil,  WT,  is  the  linear  combination  of  the  wavefront  error 


from  each  of  the  propagations  and  the  lens  transfonnation 

WT  —  Wl+A  +  W2. 


(7.15) 


Finally,  the  total  wavefront  is  decomposed  to  determine  the  each  of  the  Zernike 
coefficients,  Z.,  the  using  the  orthogonality  of  the  each  Zemike  polynomials,  4>n  such 


that 


(7.16) 
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